Context Navigation

source: trunk/source/processes/hadronic/models/cascade/cascade/src/G4Fissioner.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.2 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: G4Fissioner.cc,v 1.28 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	// 20100318 M. Kelsey -- Bug fix setting mass with G4LV
30	// 20100319 M. Kelsey -- Use new generateWithRandomAngles for theta,phi stuff;
31	// eliminate some unnecessary std::pow()
32	// 20100413 M. Kelsey -- Pass G4CollisionOutput by ref to ::collide()
33	// 20100517 M. Kelsey -- Inherit from common base class
34
35	#include "G4Fissioner.hh"
36	#include "G4CollisionOutput.hh"
37	#include "G4InuclNuclei.hh"
38	#include "G4InuclParticle.hh"
39	#include "G4FissionStore.hh"
40	#include "G4FissionConfiguration.hh"
41	#include "G4NucleiProperties.hh"
42	#include "G4HadTmpUtil.hh"
43	#include "G4InuclSpecialFunctions.hh"
44
45	using namespace G4InuclSpecialFunctions;
46
47
48	G4Fissioner::G4Fissioner() : G4VCascadeCollider("G4Fissioner") {}
49
50	void G4Fissioner::collide(G4InuclParticle* /bullet/,
51	G4InuclParticle* target,
52	G4CollisionOutput& output) {
53	if (verboseLevel > 3) {
54	G4cout << " >>> G4Fissioner::collide" << G4endl;
55	}
56
57	// const G4int itry_max = 1000;
58
59	if (G4InuclNuclei* nuclei_target = dynamic_cast<G4InuclNuclei*>(target)) {
60
61	if (verboseLevel > 3) {
62	G4cout << " Fissioner input " << G4endl;
63
64	nuclei_target->printParticle();
65	}
66
67	G4double A = nuclei_target->getA();
68	G4double Z = nuclei_target->getZ();
69	G4double EEXS = nuclei_target->getExitationEnergy();
70	G4double mass_in = nuclei_target->getMass();
71	G4double e_in = mass_in + 0.001 * EEXS;
72	G4double PARA = 0.055 * G4cbrt(AA) (G4cbrt(A - Z) + G4cbrt(Z));
73	G4double TEM = std::sqrt(EEXS / PARA);
74	G4double TETA = 0.494 * G4cbrt(A) * TEM;
75
76	TETA = TETA / std::sinh(TETA);
77
78	if (A < 246.0) PARA += (nucleiLevelDensity(A) - PARA) * TETA;
79
80	G4double A1 = G4int(A / 2.0 + 1.1);
81	G4double Z1;
82	G4double A2 = A - A1;
83	G4double ALMA = -1000.0;
84	// G4double DM1 = bindingEnergy(A, Z);
85	G4double DM1 = G4NucleiProperties::GetBindingEnergy(G4lrint(A), G4lrint(Z));
86	G4double EVV = EEXS - DM1;
87	G4double DM2 = bindingEnergyAsymptotic(A, Z);
88	G4double DTEM = (A < 220.0 ? 0.5 : 1.15);
89
90	TEM += DTEM;
91
92	std::vector<G4double> AL1(2, -0.15);
93	std::vector<G4double> BET1(2, 0.05);
94	G4FissionStore fissionStore;
95	G4double R12 = G4cbrt(A1) + G4cbrt(A2);
96
97	for (G4int i = 0; i < 50 && A1 > 30.0; i++) {
98	A1 -= 1.0;
99	A2 = A - A1;
100	G4double X3 = 1.0 / G4cbrt(A1);
101	G4double X4 = 1.0 / G4cbrt(A2);
102	Z1 = G4int(getZopt(A1, A2, Z, X3, X4, R12)) - 1.0;
103	std::vector<G4double> EDEF1(2);
104	G4double Z2 = Z - Z1;
105	G4double VPOT, VCOUL;
106
107	potentialMinimization(VPOT, EDEF1, VCOUL, A1, A2, Z1, Z2, AL1, BET1, R12);
108
109	// G4double DM3 = bindingEnergy(A1, Z1);
110	G4double DM3 = G4NucleiProperties::GetBindingEnergy(G4lrint(A1), G4lrint(Z1));
111	G4double DM4 = bindingEnergyAsymptotic(A1, Z1);
112	// G4double DM5 = bindingEnergy(A2, Z2);
113	G4double DM5 = G4NucleiProperties::GetBindingEnergy(G4lrint(A2), G4lrint(Z2));
114	G4double DM6 = bindingEnergyAsymptotic(A2, Z2);
115	G4double DMT1 = DM4 + DM6 - DM2;
116	G4double DMT = DM3 + DM5 - DM1;
117	G4double EZL = EEXS + DMT - VPOT;
118
119	if(EZL > 0.0) { // generate fluctuations
120	// faster, using randomGauss
121	G4double C1 = std::sqrt(getC2(A1, A2, X3, X4, R12) / TEM);
122	G4double DZ = randomGauss(C1);
123
124	DZ = DZ > 0.0 ? G4int(DZ + 0.5) : -G4int(std::fabs(DZ - 0.5));
125	Z1 += DZ;
126	Z2 -= DZ;
127
128	G4double DEfin = randomGauss(TEM);
129	G4double EZ = (DMT1 + (DMT - DMT1) * TETA - VPOT + DEfin) / TEM;
130
131	if (EZ >= ALMA) ALMA = EZ;
132	G4double EK = VCOUL + DEfin + 0.5 * TEM;
133	// G4double EV = EVV + bindingEnergy(A1, Z1) + bindingEnergy(A2, Z2) - EK;
134	G4double EV = EVV
135	+ G4NucleiProperties::GetBindingEnergy(G4lrint(A1), G4lrint(Z1))
136	+ G4NucleiProperties::GetBindingEnergy(G4lrint(A2), G4lrint(Z2))
137	- EK;
138
139	if (EV > 0.0) fissionStore.addConfig(A1, Z1, EZ, EK, EV);
140	};
141	};
142
143	G4int store_size = fissionStore.size();
144
145	if (store_size > 0) {
146
147	G4FissionConfiguration config =
148	fissionStore.generateConfiguration(ALMA, inuclRndm());
149
150	A1 = config.afirst;
151	A2 = A - A1;
152	Z1 = config.zfirst;
153
154	G4double Z2 = Z - Z1;
155	G4InuclNuclei nuclei1(A1, Z1);
156	nuclei1.setModel(7); // sign in the modelId (=G4Fissioner)
157	G4InuclNuclei nuclei2(A2, Z2);
158	nuclei2.setModel(7);
159
160	G4double mass1 = nuclei1.getMass();
161	G4double mass2 = nuclei2.getMass();
162	G4double EK = config.ekin;
163	G4double pmod = std::sqrt(0.001 * EK * mass1 * mass2 / mass_in);
164
165	G4LorentzVector mom1 = generateWithRandomAngles(pmod, mass1);
166	G4LorentzVector mom2; mom2.setVectM(-mom1.vect(), mass2);
167
168	G4double e_out = mom1.e() + mom2.e();
169	G4double EV = 1000.0 * (e_in - e_out) / A;
170
171	if (EV > 0.0) {
172	G4double EEXS1 = EV*A1;
173	G4double EEXS2 = EV*A2;
174	G4InuclNuclei nuclei1(mom1, A1, Z1);
175
176	nuclei1.setModel(7);
177	nuclei1.setExitationEnergy(EEXS1);
178	nuclei1.setEnergy();
179	output.addTargetFragment(nuclei1);
180
181	G4InuclNuclei nuclei2(mom2, A2, Z2);
182
183	nuclei2.setModel(7);
184	nuclei2.setExitationEnergy(EEXS2);
185	nuclei2.setEnergy();
186	output.addTargetFragment(nuclei2);
187
188	if (verboseLevel > 3) {
189	nuclei1.printParticle();
190	nuclei2.printParticle();
191	}
192	};
193	};
194
195	} else {
196	G4cout << " Fissioner -> target is not nuclei " << G4endl;
197	};
198
199	return;
200	}
201
202	G4double G4Fissioner::getC2(G4double A1,
203	G4double A2,
204	G4double X3,
205	G4double X4,
206	G4double R12) const {
207
208	if (verboseLevel > 3) {
209	G4cout << " >>> G4Fissioner::getC2" << G4endl;
210	}
211
212	G4double C2 = 124.57 * (1.0 / A1 + 1.0 / A2) + 0.78 * (X3 + X4) - 176.9 *
213	((X3X3X3X3) + (X4X4X4X4)) + 219.36 * (1.0 / (A1 * A1) + 1.0 / (A2 * A2)) - 1.108 / R12;
214
215	return C2;
216	}
217
218	G4double G4Fissioner::getZopt(G4double A1,
219	G4double A2,
220	G4double ZT,
221	G4double X3,
222	G4double X4,
223	G4double R12) const {
224
225	if (verboseLevel > 3) {
226	G4cout << " >>> G4Fissioner::getZopt" << G4endl;
227	}
228
229	G4double Zopt = (87.7 * (X4 - X3) * (1.0 - 1.25 * (X4 + X3)) +
230	ZT * ((124.57 / A2 + 0.78 * X4 - 176.9 * (X4X4X4X4) + 219.36 / (A2 A2)) - 0.554 / R12)) /
231	getC2(A1, A2, X3, X4, R12);
232
233	return Zopt;
234	}
235
236	void G4Fissioner::potentialMinimization(G4double& VP,
237	std::vector<G4double> & ED,
238	G4double& VC,
239	G4double AF,
240	G4double AS,
241	G4double ZF,
242	G4double ZS,
243	std::vector<G4double>& AL1,
244	std::vector<G4double>& BET1,
245	G4double& R12) const {
246
247	if (verboseLevel > 3) {
248	G4cout << " >>> G4Fissioner::potentialMinimization" << G4endl;
249	}
250
251	const G4double huge_num = 2.0e35;
252	const G4int itry_max = 2000;
253	const G4double DSOL1 = 1.0e-6;
254	const G4double DS1 = 0.3;
255	const G4double DS2 = 1.0 / DS1 / DS1;
256	G4double A1[2];
257	A1[0] = AF;
258	A1[1] = AS;
259	G4double Z1[2];
260	Z1[0] = ZF;
261	Z1[1] = ZS;
262	G4double D = 1.01844 * ZF * ZS;
263	G4double D0 = 1.0e-3 * D;
264	G4double R[2];
265	R12 = 0.0;
266	G4double C[2];
267	G4double F[2];
268	G4double Y1;
269	G4double Y2;
270	G4int i;
271
272	for (i = 0; i < 2; i++) {
273	R[i] = G4cbrt(A1[i]);
274	Y1 = R[i] * R[i];
275	Y2 = Z1[i] * Z1[i] / R[i];
276	C[i] = 6.8 * Y1 - 0.142 * Y2;
277	F[i] = 12.138 * Y1 - 0.145 * Y2;
278	};
279
280	G4double SAL[2];
281	G4double SBE[2];
282	G4double X[2];
283	G4double X1[2];
284	G4double X2[2];
285	G4double RAL[2];
286	G4double RBE[2];
287	G4double A[4][4];
288	G4double B[4];
289	G4int itry = 0;
290
291	while (itry < itry_max) {
292	itry++;
293	G4double S = 0.0;
294
295	for (i = 0; i < 2; i++) {
296	S += R[i] * (1.0 + AL1[i] + BET1[i] - 0.257 * AL1[i] * BET1[i]);
297	};
298	R12 = 0.0;
299	Y1 = 0.0;
300	Y2 = 0.0;
301
302	for (i = 0; i < 2; i++) {
303	SAL[i] = R[i] * (1.0-0.257 * BET1[i]);
304	SBE[i] = R[i] * (1.0-0.257 * AL1[i]);
305	X[i] = R[i] / S;
306	X1[i] = X[i] * X[i];
307	X2[i] = X[i] * X1[i];
308	Y1 += AL1[i] * X1[i];
309	Y2 += BET1[i] * X2[i];
310	R12 += R[i] * (1.0 - AL1[i] * (1.0 - 0.6 * X[i]) + BET1[i] * (1.0 - 0.429 * X1[i]));
311	};
312
313	G4double Y3 = -0.6 * Y1 + 0.857 * Y2;
314	G4double Y4 = (1.2 * Y1 - 2.571 * Y2) / S;
315	G4double R2 = D0 / (R12 * R12);
316	G4double R3 = 2.0 * R2 / R12;
317
318	for (i = 0; i < 2; i++) {
319	RAL[i] = -R[i] * (1.0 - 0.6 * X[i]) + SAL[i] * Y3;
320	RBE[i] = R[i] * (1.0 - 0.429 * X1[i]) + SBE[i] * Y3;
321	};
322
323	G4double DX1;
324	G4double DX2;
325
326	for (i = 0; i < 2; i++) {
327
328	for (G4int j = 0; j < 2; j++) {
329	G4double DEL1 = i == j ? 1.0 : 0.0;
330	DX1 = 0.0;
331	DX2 = 0.0;
332
333	if (std::fabs(AL1[i]) >= DS1) {
334	G4double XXX = AL1[i] * AL1[i] * DS2;
335	G4double DEX = XXX > 100.0 ? huge_num : std::exp(XXX);
336	DX1 = 2.0 * (1.0 + 2.0 * AL1[i] * AL1[i] * DS2) * DEX * DS2;
337	};
338
339	if (std::fabs(BET1[i]) >= DS1) {
340	G4double XXX = BET1[i] * BET1[i] * DS2;
341	G4double DEX = XXX > 100.0 ? huge_num : std::exp(XXX);
342	DX2 = 2.0 * (1.+2.0 * BET1[i] * BET1[i] * DS2) * DEX * DS2;
343	};
344
345	G4double DEL = 2.0e-3 * DEL1;
346	A[i][j] = R3 * RBE[i] * RBE[j] -
347	R2 * (-0.6 *
348	(X1[i] * SAL[j] +
349	X1[j] * SAL[i]) + SAL[i] * SAL[j] * Y4) +
350	DEL * C[i] + DEL1 * DX1;
351	G4int i1 = i + 2;
352	G4int j1 = j + 2;
353	A[i1][j1] = R3 * RBE[i] * RBE[j]
354	- R2 * (0.857 *
355	(X2[i] * SBE[j] +
356	X2[j] * SBE[i]) + SBE[i] * SBE[j] * Y4) +
357	DEL * F[i] + DEL1 * DX2;
358	A[i][j1] = R3 * RAL[i] * RBE[j] -
359	R2 * (0.857 *
360	(X2[j] * SAL[i] -
361	0.6 * X1[i] * SBE[j]) + SBE[j] * SAL[i] * Y4 -
362	0.257 * R[i] * Y3 * DEL1);
363	A[j1][i] = A[i][j1];
364	};
365	};
366
367	for (i = 0; i < 2; i++) {
368	DX1 = 0.0;
369	DX2 = 0.0;
370
371	if (std::fabs(AL1[i]) >= DS1) DX1 = 2.0 * AL1[i] * DS2 * std::exp(AL1[i] * AL1[i] * DS2);
372
373	if (std::fabs(BET1[i]) >= DS1) DX2 = 2.0 * BET1[i] * DS2 * std::exp(BET1[i] * BET1[i] * DS2);
374	B[i] = R2 * RAL[i] - 2.0e-3 * C[i] * AL1[i] + DX1;
375	B[i + 2] = R2 * RBE[i] - 2.0e-3 * F[i] * BET1[i] + DX2;
376	};
377
378	G4double ST = 0.0;
379	G4double ST1 = 0.0;
380
381	for (i = 0; i < 4; i++) {
382	ST += B[i] * B[i];
383
384	for (G4int j = 0; j < 4; j++) ST1 += A[i][j] * B[i] * B[j];
385	};
386
387	G4double STEP = ST / ST1;
388	G4double DSOL = 0.0;
389
390	for (i = 0; i < 2; i++) {
391	AL1[i] += B[i] * STEP;
392	BET1[i] += B[i + 2] * STEP;
393	DSOL += B[i] * B[i] + B[i + 2] * B[i + 2];
394	};
395	DSOL = std::sqrt(DSOL);
396
397	if (DSOL < DSOL1) break;
398	};
399
400	if (verboseLevel > 3) {
401	if (itry == itry_max)
402	G4cout << " maximal number of iterations in potentialMinimization " << G4endl
403	<< " A1 " << AF << " Z1 " << ZF << G4endl;
404
405	};
406
407	for (i = 0; i < 2; i++) ED[i] = F[i] * BET1[i] * BET1[i] + C[i] * AL1[i] * AL1[i];
408
409	VC = D / R12;
410	VP = VC + ED[0] + ED[1];
411	}

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

Download in other formats: