Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

G4Abla.cc @ 1350

Last change on this file since 1350 was 1350, checked in by garnier, 14 years ago
update to last version 4.9.4
File size: 188.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	//
26	// $Id: G4Abla.cc,v 1.1 2008/02/27 18:31:11 miheikki Exp $
27	// Translation of INCL4.2/ABLA V3
28	// Pekka Kaitaniemi, HIP (translation)
29	// Christelle Schmidt, IPNL (fission code)
30	// Alain Boudard, CEA (contact person INCL/ABLA)
31	// Aatos Heikkinen, HIP (project coordination)
32
33	#include <time.h>
34
35	#include "G4Abla.hh"
36	#include "G4InclAblaDataFile.hh"
37	#include "Randomize.hh"
38	#include <assert.h>
39
40	G4Abla::G4Abla()
41	{
42	ilast = 0;
43	}
44
45	G4Abla::G4Abla(G4Hazard hazard, G4Volant volant)
46	{
47	verboseLevel = 0;
48	ilast = 0;
49	volant = volant; // ABLA internal particle data
50	volant->iv = 0;
51	hazard = hazard; // Random seeds
52
53	varntp = new G4VarNtp();
54	pace = new G4Pace();
55	ald = new G4Ald();
56	ablamain = new G4Ablamain();
57	emdpar = new G4Emdpar();
58	eenuc = new G4Eenuc();
59	ec2sub = new G4Ec2sub();
60	ecld = new G4Ecld();
61	fb = new G4Fb();
62	fiss = new G4Fiss();
63	opt = new G4Opt();
64	}
65
66	G4Abla::G4Abla(G4Hazard aHazard, G4Volant aVolant, G4VarNtp *aVarntp)
67	{
68	verboseLevel = 0;
69	ilast = 0;
70	volant = aVolant; // ABLA internal particle data
71	volant->iv = 0;
72	hazard = aHazard; // Random seeds
73	varntp = aVarntp; // Output data structure
74	varntp->ntrack = 0;
75
76	pace = new G4Pace();
77	ald = new G4Ald();
78	ablamain = new G4Ablamain();
79	emdpar = new G4Emdpar();
80	eenuc = new G4Eenuc();
81	ec2sub = new G4Ec2sub();
82	ecld = new G4Ecld();
83	fb = new G4Fb();
84	fiss = new G4Fiss();
85	opt = new G4Opt();
86	}
87
88	G4Abla::~G4Abla()
89	{
90	delete pace;
91	delete ald;
92	delete ablamain;
93	delete emdpar;
94	delete eenuc;
95	delete ec2sub;
96	delete ecld;
97	delete fb;
98	delete fiss;
99	delete opt;
100	}
101
102	// Main interface to the evaporation
103
104	// Possible problem with generic Geant4 interface: ABLA evaporation
105	// needs angular momentum information (calculated by INCL) to
106	// work. Maybe there is a way to obtain this information from
107	// G4Fragment?
108
109	void G4Abla::breakItUp(G4double nucleusA, G4double nucleusZ, G4double nucleusMass, G4double excitationEnergy,
110	G4double angularMomentum, G4double recoilEnergy, G4double momX, G4double momY, G4double momZ,
111	G4int eventnumber)
112	{
113	const G4double uma = 931.4942;
114	const G4double melec = 0.511;
115	const G4double fmp = 938.27231;
116	const G4double fmn = 939.56563;
117
118	G4double alrem = 0.0, berem = 0.0, garem = 0.0;
119	G4double R[4][4]; // Rotation matrix
120	G4double csdir1[4];
121	G4double csdir2[4];
122	G4double csrem[4];
123	G4double pfis_rem[4];
124	G4double pf1_rem[4];
125	for(G4int init_i = 0; init_i < 4; init_i++) {
126	csdir1[init_i] = 0.0;
127	csdir2[init_i] = 0.0;
128	csrem[init_i] = 0.0;
129	pfis_rem[init_i] = 0.0;
130	pf1_rem[init_i] = 0.0;
131	for(G4int init_j = 0; init_j < 4; init_j++) {
132	R[init_i][init_j] = 0.0;
133	}
134	}
135
136	G4double plab1 = 0.0, gam1 = 0.0, eta1 = 0.0;
137	G4double plab2 = 0.0, gam2 = 0.0, eta2 = 0.0;
138
139	G4double sitet = 0.0;
140	G4double stet1 = 0.0;
141	G4double stet2 = 0.0;
142
143	G4int nbpevap = 0;
144	G4int mempaw = 0, memiv = 0;
145
146	G4double e_evapo = 0.0;
147	G4double el = 0.0;
148	G4double fmcv = 0.0;
149
150	G4double aff1 = 0.0;
151	G4double zff1 = 0.0;
152	G4double eff1 = 0.0;
153	G4double aff2 = 0.0;
154	G4double zff2 = 0.0;
155	G4double eff2 = 0.0;
156
157	G4double v1 = 0.0, v2 = 0.0;
158
159	G4double t2 = 0.0;
160	G4double ctet1 = 0.0;
161	G4double ctet2 = 0.0;
162	G4double phi1 = 0.0;
163	G4double phi2 = 0.0;
164	G4double p2 = 0.0;
165	G4double epf2_out = 0.0 ;
166	G4int lma_pf1 = 0, lmi_pf1 = 0;
167	G4int lma_pf2 = 0, lmi_pf2 = 0;
168	G4int nopart = 0;
169
170	G4double cst = 0.0, sst = 0.0, csf = 0.0, ssf = 0.0;
171
172	G4double zf = 0.0, af = 0.0, mtota = 0.0, pleva = 0.0, pxeva = 0.0, pyeva = 0.0;
173	G4int ff = 0;
174	G4int inum = eventnumber;
175	G4int inttype = 0;
176	G4double esrem = excitationEnergy;
177
178	G4double aprf = nucleusA;
179	G4double zprf = nucleusZ;
180	G4double mcorem = nucleusMass;
181	G4double ee = excitationEnergy;
182	G4double jprf = angularMomentum; // actually root-mean-squared
183
184	G4double erecrem = recoilEnergy;
185	G4double trem = 0.0;
186	G4double pxrem = momX;
187	G4double pyrem = momY;
188	G4double pzrem = momZ;
189
190	G4double remmass = 0.0;
191
192	varntp->ntrack = 0;
193	// volant->iv = 0;
194	volant->iv = 1;
195
196	G4double pcorem = std::sqrt(erecrem(erecrem +2.938.2796*nucleusA));
197	// G4double pcorem = std::sqrt(std::pow(momX,2) + std::pow(momY,2) + std::pow(momZ,2));
198	// assert(isnan(pcorem) == false);
199	if(esrem >= 1.0e-3) {
200	evapora(zprf,aprf,ee,jprf, &zf, &af, &mtota, &pleva, &pxeva, &pyeva, &ff, &inttype, &inum);
201	}
202	else {
203	ff = 0;
204	zf = zprf;
205	af = aprf;
206	pxeva = pxrem;
207	pyeva = pyrem;
208	pleva = pzrem;
209	}
210	// assert(isnan(zf) == false);
211	// assert(isnan(af) == false);
212	// assert(isnan(ee) == false);
213
214	if (ff == 1) {
215	// Fission:
216	// variable ee: Energy of fissioning nucleus above the fission barrier.
217	// Calcul des impulsions des particules evaporees (avant fission)
218	// dans le systeme labo.
219
220	trem = double(erecrem);
221	remmass = pace2(aprf,zprf) + aprfuma - zprfmelec; // canonic
222	// remmass = mcorem + double(esrem); // ok
223	// remmass = mcorem; //cugnon
224	varntp->kfis = 1;
225	G4double gamrem = (remmass + trem)/remmass;
226	G4double etrem = std::sqrt(trem(trem + 2.0remmass))/remmass;
227	// assert(isnan(etrem) == false);
228
229	// This is not treated as accurately as for the non fission case for which
230	// the remnant mass is computed to satisfy the energy conservation
231	// of evaporated particles. But it is not bad and more canonical!
232	remmass = pace2(aprf,zprf) + aprfuma - zprfmelec+double(esrem); // !canonic
233	// Essais avec la masse de KHS (9/2002):
234	el = 0.0;
235	mglms(aprf,zprf,0,&el);
236	remmass = zprffmp + (aprf-zprf)fmn + el + double(esrem);
237
238	gamrem = std::sqrt(std::pow(pcorem,2) + std::pow(remmass,2))/remmass;
239	// assert(isnan(gamrem) == false);
240	etrem = pcorem/remmass;
241	// assert(isnan(etrem) == false);
242
243	alrem = pxrem/pcorem;
244	// assert(isnan(alrem) == false);
245	berem = pyrem/pcorem;
246	// assert(isnan(berem) == false);
247	garem = pzrem/pcorem;
248	// assert(isnan(garem) == false);
249
250	csrem[0] = 0.0; // Should not be used.
251	csrem[1] = alrem;
252	csrem[2] = berem;
253	csrem[3] = garem;
254
255	// C Pour VÃ©rif Remnant = evapo(Pre fission) + Noyau_fissionant (systÃšme Remnant)
256	G4double bil_e = 0.0;
257	G4double bil_px = 0.0;
258	G4double bil_py = 0.0;
259	G4double bil_pz = 0.0;
260	G4double masse = 0.0;
261
262	for(G4int iloc = 1; iloc <= volant->iv; iloc++) { //DO iloc=1,iv
263	// assert(isnan(volant->zpcv[iloc]) == false);
264	// assert(volant->acv[iloc] != 0);
265	// assert(volant->zpcv[iloc] != 0);
266	mglms(double(volant->acv[iloc]),double(volant->zpcv[iloc]),0,&el);
267	// assert(isnan(el) == false);
268	masse = volant->zpcv[iloc]fmp + (volant->acv[iloc] - volant->zpcv[iloc])fmn + el;
269	// assert(isnan(masse) == false);
270	bil_e = bil_e + std::sqrt(std::pow(volant->pcv[iloc],2) + std::pow(masse,2));
271	// assert(isnan(bil_e) == false);
272	bil_px = bil_px + volant->pcv[iloc]*(volant->xcv[iloc]);
273	bil_py = bil_py + volant->pcv[iloc]*(volant->ycv[iloc]);
274	bil_pz = bil_pz + volant->pcv[iloc]*(volant->zcv[iloc]);
275	} // enddo
276	// C Ce bilan (impulsion nulle) est parfait. (Bil_Px=Bil_Px+PXEVA....)
277
278	G4int ndec = 1;
279
280	if(volant->iv != 0) { //then
281	if(verboseLevel > 2) {
282	G4cout <<"varntp->ntrack = " << varntp->ntrack << G4endl;
283	G4cout <<"1st Translab: Adding indices from " << ndec << " to " << volant->iv << G4endl;
284	}
285	nopart = varntp->ntrack - 1;
286	translab(gamrem,etrem,csrem,nopart,ndec);
287	if(verboseLevel > 2) {
288	G4cout <<"Translab complete!" << G4endl;
289	G4cout <<"varntp->ntrack = " << varntp->ntrack << G4endl;
290	}
291	}
292	nbpevap = volant->iv; // nombre de particules d'evaporation traitees
293
294	// C
295	// C Now calculation of the fission fragment distribution including
296	// C evaporation from the fragments.
297	// C
298
299	// C Distribution of the fission fragments:
300
301	// void fissionDistri(G4double a,G4double z,G4double e,
302	// G4double &a1,G4double &z1,G4double &e1,G4double &v1,
303	// G4double &a2,G4double &z2,G4double &e2,G4double &v2);
304
305	fissionDistri(af,zf,ee,aff1,zff1,eff1,v1,aff2,zff2,eff2,v2);
306
307	// C verif des A et Z decimaux:
308	G4int na_f = int(std::floor(af + 0.5));
309	G4int nz_f = int(std::floor(zf + 0.5));
310	varntp->izfis = nz_f; // copie dans le ntuple
311	varntp->iafis = na_f;
312	G4int na_pf1 = int(std::floor(aff1 + 0.5));
313	G4int nz_pf1 = int(std::floor(zff1 + 0.5));
314	G4int na_pf2 = int(std::floor(aff2 + 0.5));
315	G4int nz_pf2 = int(std::floor(zff2 + 0.5));
316
317	if((na_f != (na_pf1+na_pf2)) \|\| (nz_f != (nz_pf1+nz_pf2))) {
318	if(verboseLevel > 2) {
319	G4cout <<"problemes arrondis dans la fission " << G4endl;
320	G4cout << "af,zf,aff1,zff1,aff2,zff2" << G4endl;
321	G4cout << af <<" , " << zf <<" , " << aff1 <<" , " << zff1 <<" , " << aff2 <<" , " << zff2 << G4endl;
322	G4cout << "a,z,a1,z1,a2,z2 integer" << G4endl;
323	G4cout << na_f <<" , " << nz_f <<" , " << na_pf1 <<" , " << nz_pf1 <<" , " << na_pf2 <<" , " << nz_pf2 << G4endl;
324	}
325	}
326
327	// Calcul de l'impulsion des PF dans le syteme noyau de fission:
328	G4int kboud = idnint(zf);
329	G4int jboud = idnint(af-zf);
330	//G4double ef = fb->efa[kboud][jboud]; // barriere de fission
331	G4double ef = fb->efa[jboud][kboud]; // barriere de fission
332	// assert(isnan(ef) == false);
333	varntp->estfis = ee + ef; // copie dans le ntuple
334
335	// C MASSEF = pace2(AF,ZF)
336	// C MASSEF = MASSEF + AFUMA - ZFMELEC + EE + EF
337	// C MASSE1 = pace2(DBLE(AFF1),DBLE(ZFF1))
338	// C MASSE1 = MASSE1 + AFF1UMA - ZFF1MELEC + EFF1
339	// C MASSE2 = pace2(DBLE(AFF2),DBLE(ZFF2))
340	// C MASSE2 = MASSE2 + AFF2UMA - ZFF2MELEC + EFF2
341	// C WRITE(6,*) 'MASSEF,MASSE1,MASSE2',MASSEF,MASSE1,MASSE2
342	// C MGLMS est la fonction de masse cohÃ©rente avec KHS evapo-fis.
343	// C Attention aux parametres, ici 0=OPTSHP, NO microscopic correct.
344	mglms(af,zf,0,&el);
345	// assert(isnan(el) == false);
346	G4double massef = zffmp + (af - zf)fmn + el + ee + ef;
347	// assert(isnan(massef) == false);
348	mglms(double(aff1),double(zff1),0,&el);
349	// assert(isnan(el) == false);
350	G4double masse1 = zff1fmp + (aff1-zff1)fmn + el + eff1;
351	// assert(isnan(masse1) == false);
352	mglms(aff2,zff2,0,&el);
353	// assert(isnan(el) == false);
354	G4double masse2 = zff2fmp + (aff2 - zff2)fmn + el + eff2;
355	// assert(isnan(masse2) == false);
356	// C WRITE(6,*) 'MASSEF,MASSE1,MASSE2',MASSEF,MASSE1,MASSE2
357	G4double b = massef - masse1 - masse2;
358	if(b < 0.0) { //then
359	b=0.0;
360	if(verboseLevel > 2) {
361	G4cout <<"anomalie dans la fission: " << G4endl;
362	G4cout << inum<< " , " << af<< " , " <<zf<< " , " <<massef<< " , " <<aff1<< " , " <<zff1<< " , " <<masse1<< " , " <<aff2<< " , " <<zff2<< " , " << masse2 << G4endl;
363	}
364	} //endif
365	G4double t1 = b(b + 2.0masse2)/(2.0*massef);
366	// assert(isnan(t1) == false);
367	G4double p1 = std::sqrt(t1(t1 + 2.0masse1));
368	// assert(isnan(p1) == false);
369
370	G4double rndm;
371	standardRandom(&rndm, &(hazard->igraine[13]));
372	ctet1 = 2.0*rndm - 1.0;
373	standardRandom(&rndm,&(hazard->igraine[9]));
374	phi1 = rndm2.03.141592654;
375
376	// C ----Coefs de la transformation de Lorentz (noyau de fission -> Remnant)
377	G4double peva = std::pow(pxeva,2) + std::pow(pyeva,2) + std::pow(pleva,2);
378	G4double gamfis = std::sqrt(std::pow(massef,2) + peva)/massef;
379	// assert(isnan(gamfis) == false);
380	peva = std::sqrt(peva);
381	// assert(isnan(peva) == false);
382	G4double etfis = peva/massef;
383
384	G4double epf1_in = 0.0;
385	G4double epf1_out = 0.0;
386
387	// C ----Matrice de rotation (noyau de fission -> Remnant)
388	if(peva >= 1.0e-4) {
389	sitet = std::sqrt(std::pow(pxeva,2)+std::pow(pyeva,2))/peva;
390	// assert(isnan(sitet) == false);
391	}
392	if(sitet > 1.0e-5) { //then
393	G4double cstet = pleva/peva;
394	G4double siphi = pyeva/(sitet*peva);
395	G4double csphi = pxeva/(sitet*peva);
396
397	R[1][1] = cstet*csphi;
398	R[1][2] = -siphi;
399	R[1][3] = sitet*csphi;
400	R[2][1] = cstet*siphi;
401	R[2][2] = csphi;
402	R[2][3] = sitet*siphi;
403	R[3][1] = -sitet;
404	R[3][2] = 0.0;
405	R[3][3] = cstet;
406	}
407	else {
408	R[1][1] = 1.0;
409	R[1][2] = 0.0;
410	R[1][3] = 0.0;
411	R[2][1] = 0.0;
412	R[2][2] = 1.0;
413	R[2][3] = 0.0;
414	R[3][1] = 0.0;
415	R[3][2] = 0.0;
416	R[3][3] = 1.0;
417	} // endif
418	// c test de verif:
419
420	if((zff1 <= 0.0) \|\| (aff1 <= 0.0) \|\| (aff1 < zff1)) { //then
421	if(verboseLevel > 2) {
422	G4cout <<"zf = " << zf <<" af = " << af <<"ee = " << ee <<"zff1 = " << zff1 <<"aff1 = " << aff1 << G4endl;
423	}
424	}
425	else {
426	// C ---------------------- PF1 will evaporate
427	epf1_in = double(eff1);
428	epf1_out = epf1_in;
429	// void evapora(G4double zprf, G4double aprf, G4double ee, G4double jprf,
430	// G4double zf_par, G4double af_par, G4double *mtota_par,
431	// G4double pleva_par, G4double pxeva_par, G4double *pyeva_par,
432	// G4double ff_par, G4int inttype_par, G4int *inum_par);
433	G4double zf1, af1, malpha1, ffpleva1, ffpxeva1, ffpyeva1;
434	G4int ff1, ftype1;
435	evapora(zff1, aff1, epf1_out, 0.0, &zf1, &af1, &malpha1, &ffpleva1,
436	&ffpxeva1, &ffpyeva1, &ff1, &ftype1, &inum);
437	// C On ajoute le fragment:
438	// assert(af1 > 0);
439	volant->iv = volant->iv + 1;
440	// assert(af1 != 0);
441	// assert(zf1 != 0);
442	volant->acv[volant->iv] = af1;
443	volant->zpcv[volant->iv] = zf1;
444	if(verboseLevel > 2) {
445	G4cout <<"Added fission fragment: a = " << volant->acv[volant->iv] << " z = " << volant->zpcv[volant->iv] << G4endl;
446	}
447	peva = std::sqrt(std::pow(ffpxeva1,2) + std::pow(ffpyeva1,2) + std::pow(ffpleva1,2));
448	// assert(isnan(peva) == false);
449	volant->pcv[volant->iv] = peva;
450	if(peva > 0.001) { // then
451	volant->xcv[volant->iv] = ffpxeva1/peva;
452	volant->ycv[volant->iv] = ffpyeva1/peva;
453	volant->zcv[volant->iv] = ffpleva1/peva;
454	}
455	else {
456	volant->xcv[volant->iv] = 1.0;
457	volant->ycv[volant->iv] = 0.0;
458	volant->zcv[volant->iv] = 0.0;
459	} // end if
460
461	// C Pour VÃ©rif evapo de PF1 dans le systeme du Noyau Fissionant
462	G4double bil1_e = 0.0;
463	G4double bil1_px = 0.0;
464	G4double bil1_py=0.0;
465	G4double bil1_pz=0.0;
466	for(G4int iloc = nbpevap + 1; iloc <= volant->iv; iloc++) { //do iloc=nbpevap+1,iv
467	// for(G4int iloc = nbpevap + 1; iloc <= volant->iv + 1; iloc++) { //do iloc=nbpevap+1,iv
468	mglms(volant->acv[iloc], volant->zpcv[iloc],0,&el);
469	masse = volant->zpcv[iloc]fmp + (volant->acv[iloc] - volant->zpcv[iloc])fmn + el;
470	// assert(isnan(masse) == false);
471	bil1_e = bil1_e + std::sqrt(std::pow(volant->pcv[iloc],2) + std::pow(masse,2));
472	// assert(isnan(bil1_e) == false);
473	bil1_px = bil1_px + volant->pcv[iloc]*(volant->xcv[iloc]);
474	bil1_py = bil1_py + volant->pcv[iloc]*(volant->ycv[iloc]);
475	bil1_pz = bil1_pz + volant->pcv[iloc]*(volant->zcv[iloc]);
476	} // enddo
477
478	// Calcul des cosinus directeurs de PF1 dans le Remnant et calcul
479	// des coefs pour la transformation de Lorentz Systeme PF --> Systeme Remnant
480	translabpf(masse1,t1,p1,ctet1,phi1,gamfis,etfis,R,&plab1,&gam1,&eta1,csdir1);
481
482	// calcul des impulsions des particules evaporees dans le systeme Remnant:
483	if(verboseLevel > 2) {
484	G4cout <<"2nd Translab (pf1 evap): Adding indices from " << nbpevap+1 << " to " << volant->iv << G4endl;
485	}
486	nopart = varntp->ntrack - 1;
487	translab(gam1,eta1,csdir1,nopart,nbpevap+1);
488	if(verboseLevel > 2) {
489	G4cout <<"After translab call... varntp->ntrack = " << varntp->ntrack << G4endl;
490	}
491	memiv = nbpevap + 1; // memoires pour la future transformation
492	mempaw = nopart; // remnant->labo pour pf1 et pf2.
493	lmi_pf1 = nopart + nbpevap + 1; // indices min et max dans /var_ntp/
494	lma_pf1 = nopart + volant->iv; // des particules issues de pf1
495	nbpevap = volant->iv; // nombre de particules d'evaporation traitees
496	} // end if
497	// C --------------------- End of PF1 calculation
498
499	// c test de verif:
500	if((zff2 <= 0.0) \|\| (aff2 <= 0.0) \|\| (aff2 <= zff2)) { //then
501	if(verboseLevel > 2) {
502	G4cout << zf << " " << af << " " << ee << " " << zff2 << " " << aff2 << G4endl;
503	}
504	}
505	else {
506	// C ---------------------- PF2 will evaporate
507	G4double epf2_in = double(eff2);
508	G4double epf2_out = epf2_in;
509	// void evapora(G4double zprf, G4double aprf, G4double ee, G4double jprf,
510	// G4double zf_par, G4double af_par, G4double *mtota_par,
511	// G4double pleva_par, G4double pxeva_par, G4double *pyeva_par,
512	// G4double ff_par, G4int inttype_par, G4int *inum_par);
513	G4double zf2, af2, malpha2, ffpleva2, ffpxeva2, ffpyeva2;
514	G4int ff2, ftype2;
515	evapora(zff2,aff2,epf2_out,0.0,&zf2,&af2,&malpha2,&ffpleva2,
516	&ffpxeva2,&ffpyeva2,&ff2,&ftype2,&inum);
517	// C On ajoute le fragment:
518	volant->iv = volant->iv + 1;
519	volant->acv[volant->iv] = af2;
520	volant->zpcv[volant->iv] = zf2;
521	if(verboseLevel > 2) {
522	G4cout <<"Added fission fragment: a = " << volant->acv[volant->iv] << " z = " << volant->zpcv[volant->iv] << G4endl;
523	}
524	peva = std::sqrt(std::pow(ffpxeva2,2) + std::pow(ffpyeva2,2) + std::pow(ffpleva2,2));
525	// assert(isnan(peva) == false);
526	volant->pcv[volant->iv] = peva;
527	// exit(0);
528	if(peva > 0.001) { //then
529	volant->xcv[volant->iv] = ffpxeva2/peva;
530	volant->ycv[volant->iv] = ffpyeva2/peva;
531	volant->zcv[volant->iv] = ffpleva2/peva;
532	}
533	else {
534	volant->xcv[volant->iv] = 1.0;
535	volant->ycv[volant->iv] = 0.0;
536	volant->zcv[volant->iv] = 0.0;
537	} //end if
538	// C Pour VÃ©rif evapo de PF1 dans le systeme du Noyau Fissionant
539	G4double bil2_e = 0.0;
540	G4double bil2_px = 0.0;
541	G4double bil2_py = 0.0;
542	G4double bil2_pz = 0.0;
543	// for(G4int iloc = nbpevap + 1; iloc <= volant->iv; iloc++) { //do iloc=nbpevap+1,iv
544	for(G4int iloc = nbpevap + 1; iloc <= volant->iv; iloc++) { //do iloc=nbpevap+1,iv
545	mglms(volant->acv[iloc],volant->zpcv[iloc],0,&el);
546	masse = volant->zpcv[iloc]fmp + (volant->acv[iloc] - volant->zpcv[iloc])fmn + el;
547	bil2_e = bil2_e + std::sqrt(std::pow(volant->pcv[iloc],2) + std::pow(masse,2));
548	// assert(isnan(bil2_e) == false);
549	bil2_px = bil2_px + volant->pcv[iloc]*(volant->xcv[iloc]);
550	bil2_py = bil2_py + volant->pcv[iloc]*(volant->ycv[iloc]);
551	bil2_pz = bil2_pz + volant->pcv[iloc]*(volant->zcv[iloc]);
552	} //enddo
553
554	// C ----Calcul des cosinus directeurs de PF2 dans le Remnant et calcul
555	// c des coefs pour la transformation de Lorentz Systeme PF --> Systeme Remnant
556	G4double t2 = b - t1;
557	// G4double ctet2 = -ctet1;
558	ctet2 = -1.0*ctet1;
559	assert(std::fabs(ctet2) <= 1.0);
560	// assert(isnan(ctet2) == false);
561	phi2 = dmod(phi1+3.141592654,6.283185308);
562	// assert(isnan(phi2) == false);
563	G4double p2 = std::sqrt(t2(t2+2.0masse2));
564	// assert(isnan(p2) == false);
565
566	// void translabpf(G4double masse1, G4double t1, G4double p1, G4double ctet1,
567	// G4double phi1, G4double gamrem, G4double etrem, G4double R[][4],
568	// G4double plab1, G4double gam1, G4double *eta1, G4double csdir[]);
569	translabpf(masse2,t2,p2,ctet2,phi2,gamfis,etfis,R,&plab2,&gam2,&eta2,csdir2);
570	// C
571	// C calcul des impulsions des particules evaporees dans le systeme Remnant:
572	// c
573	if(verboseLevel > 2) {
574	G4cout <<"3rd Translab (pf2 evap): Adding indices from " << nbpevap+1 << " to " << volant->iv << G4endl;
575	}
576	nopart = varntp->ntrack - 1;
577	translab(gam2,eta2,csdir2,nopart,nbpevap+1);
578	lmi_pf2 = nopart + nbpevap + 1; // indices min et max dans /var_ntp/
579	lma_pf2 = nopart + volant->iv; // des particules issues de pf2
580	} // end if
581	// C --------------------- End of PF2 calculation
582
583	// C Pour vÃ©rifications: calculs du noyau fissionant et des PF dans
584	// C le systeme du remnant.
585	for(G4int iloc = 1; iloc <= 3; iloc++) { // do iloc=1,3
586	pfis_rem[iloc] = 0.0;
587	} // enddo
588	G4double efis_rem, pfis_trav[4];
589	lorab(gamfis,etfis,massef,pfis_rem,&efis_rem,pfis_trav);
590	rotab(R,pfis_trav,pfis_rem);
591
592	stet1 = std::sqrt(1.0 - std::pow(ctet1,2));
593	// assert(isnan(stet1) == false);
594	pf1_rem[1] = p1stet1std::cos(phi1);
595	pf1_rem[2] = p1stet1std::sin(phi1);
596	pf1_rem[3] = p1*ctet1;
597	G4double e1_rem;
598	lorab(gamfis,etfis,masse1+t1,pf1_rem,&e1_rem,pfis_trav);
599	rotab(R,pfis_trav,pf1_rem);
600
601	stet2 = std::sqrt(1.0 - std::pow(ctet2,2));
602	assert(std::pow(ctet2,2) >= 0.0);
603	assert(std::pow(ctet2,2) <= 1.0);
604	// assert(isnan(stet2) == false);
605
606	G4double pf2_rem[4];
607	G4double e2_rem;
608	pf2_rem[1] = p2stet2std::cos(phi2);
609	pf2_rem[2] = p2stet2std::sin(phi2);
610	pf2_rem[3] = p2*ctet2;
611	lorab(gamfis,etfis,masse2+t2,pf2_rem,&e2_rem,pfis_trav);
612	rotab(R,pfis_trav,pf2_rem);
613	// C Verif 0: Remnant = evapo_pre_fission + Noyau Fissionant
614	bil_e = remmass - efis_rem - bil_e;
615	bil_px = bil_px + pfis_rem[1];
616	bil_py = bil_py + pfis_rem[2];
617	bil_pz = bil_pz + pfis_rem[3];
618	// C Verif 1: noyau fissionant = PF1 + PF2 dans le systeme remnant
619	// G4double bilan_e = efis_rem - e1_rem - e2_rem;
620	// G4double bilan_px = pfis_rem[1] - pf1_rem[1] - pf2_rem[1];
621	// G4double bilan_py = pfis_rem[2] - pf1_rem[2] - pf2_rem[2];
622	// G4double bilan_pz = pfis_rem[3] - pf1_rem[3] - pf2_rem[3];
623	// C Verif 2: PF1 et PF2 egaux a toutes leurs particules evaporees
624	// C (Systeme remnant)
625	if((lma_pf1-lmi_pf1) != 0) { //then
626	G4double bil_e_pf1 = e1_rem - epf1_out;
627	G4double bil_px_pf1 = pf1_rem[1];
628	G4double bil_py_pf1 = pf1_rem[2];
629	G4double bil_pz_pf1 = pf1_rem[3];
630	for(G4int ipf1 = lmi_pf1; ipf1 <= lma_pf1; ipf1++) { //do ipf1=lmi_pf1,lma_pf1
631	bil_e_pf1 = bil_e_pf1 - (std::pow(varntp->plab[ipf1],2) + std::pow(varntp->enerj[ipf1],2))/(2.0*(varntp->enerj[ipf1]));
632	cst = std::cos(varntp->tetlab[ipf1]/57.2957795);
633	sst = std::sin(varntp->tetlab[ipf1]/57.2957795);
634	csf = std::cos(varntp->philab[ipf1]/57.2957795);
635	ssf = std::sin(varntp->philab[ipf1]/57.2957795);
636	bil_px_pf1 = bil_px_pf1 - varntp->plab[ipf1]sstcsf;
637	bil_py_pf1 = bil_py_pf1 - varntp->plab[ipf1]sstssf;
638	bil_pz_pf1 = bil_pz_pf1 - varntp->plab[ipf1]*cst;
639	} // enddo
640	} //endif
641
642	if((lma_pf2-lmi_pf2) != 0) { //then
643	G4double bil_e_pf2 = e2_rem - epf2_out;
644	G4double bil_px_pf2 = pf2_rem[1];
645	G4double bil_py_pf2 = pf2_rem[2];
646	G4double bil_pz_pf2 = pf2_rem[3];
647	for(G4int ipf2 = lmi_pf2; ipf2 <= lma_pf2; ipf2++) { //do ipf2=lmi_pf2,lma_pf2
648	bil_e_pf2 = bil_e_pf2 - (std::pow(varntp->plab[ipf2],2) + std::pow(varntp->enerj[ipf2],2))/(2.0*(varntp->enerj[ipf2]));
649	G4double cst = std::cos(varntp->tetlab[ipf2]/57.2957795);
650	G4double sst = std::sin(varntp->tetlab[ipf2]/57.2957795);
651	G4double csf = std::cos(varntp->philab[ipf2]/57.2957795);
652	G4double ssf = std::sin(varntp->philab[ipf2]/57.2957795);
653	bil_px_pf2 = bil_px_pf2 - varntp->plab[ipf2]sstcsf;
654	bil_py_pf2 = bil_py_pf2 - varntp->plab[ipf2]sstssf;
655	bil_pz_pf2 = bil_pz_pf2 - varntp->plab[ipf2]*cst;
656	} // enddo
657	} //endif
658	// C
659	// C ---- Transformation systeme Remnant -> systeme labo. (evapo de PF1 ET PF2)
660	// C
661	// G4double mempaw, memiv;
662	if(verboseLevel > 2) {
663	G4cout <<"4th Translab: Adding indices from " << memiv << " to " << volant->iv << G4endl;
664	}
665	translab(gamrem,etrem,csrem,mempaw,memiv);
666	// C ***************** END of fission calculations **********************
667	}
668	else {
669	// C ********************** Evapo sans fission ***************************
670	// C Here, FF=0, --> Evapo sans fission, on ajoute le fragment:
671	// C *************************************************************************
672	varntp->kfis = 0;
673	if(verboseLevel > 2) {
674	G4cout <<"Evaporation without fission" << G4endl;
675	}
676	volant->iv = volant->iv + 1;
677	volant->acv[volant->iv] = af;
678	volant->zpcv[volant->iv] = zf;
679	G4double peva = std::sqrt(std::pow(pxeva,2)+std::pow(pyeva,2)+std::pow(pleva,2));
680	// assert(isnan(peva) == false);
681	volant->pcv[volant->iv] = peva;
682	if(peva > 0.001) { //then
683	volant->xcv[volant->iv] = pxeva/peva;
684	volant->ycv[volant->iv] = pyeva/peva;
685	volant->zcv[volant->iv] = pleva/peva;
686	}
687	else {
688	volant->xcv[volant->iv] = 1.0;
689	volant->ycv[volant->iv] = 0.0;
690	volant->zcv[volant->iv] = 0.0;
691	} // end if
692
693	// C
694	// C calcul des impulsions des particules evaporees dans le systeme labo:
695	// c
696	trem = double(erecrem);
697	// C REMMASS = pace2(APRF,ZPRF) + APRFUMA - ZPRFMELEC !Canonic
698	// C REMMASS = MCOREM + DBLE(ESREM) !OK
699	remmass = mcorem; //Cugnon
700	// C GAMREM = (REMMASS + TREM)/REMMASS !OK
701	// C ETREM = DSQRT(TREM(TREM + 2.REMMASS))/REMMASS !OK
702	csrem[0] = 0.0; // Should not be used.
703	csrem[1] = alrem;
704	csrem[2] = berem;
705	csrem[3] = garem;
706
707	// for(G4int j = 1; j <= volant->iv; j++) { //do j=1,iv
708	for(G4int j = 1; j <= volant->iv; j++) { //do j=1,iv
709	if(volant->acv[j] == 0) {
710	if(verboseLevel > 2) {
711	G4cout <<"volant->acv[" << j << "] = 0" << G4endl;
712	G4cout <<"volant->iv = " << volant->iv << G4endl;
713	}
714	}
715	if(volant->acv[j] > 0) {
716	assert(volant->acv[j] != 0);
717	// assert(volant->zpcv[j] != 0);
718	mglms(volant->acv[j],volant->zpcv[j],0,&el);
719	fmcv = volant->zpcv[j]fmp + (volant->acv[j] - volant->zpcv[j])fmn + el;
720	e_evapo = e_evapo + std::sqrt(std::pow(volant->pcv[j],2) + std::pow(fmcv,2));
721	// assert(isnan(e_evapo) == false);
722	}
723	} // enddo
724
725	// C Redefinition pour conservation d'impulsion!!!
726	// C this mass obtained by energy balance is very close to the
727	// C mass of the remnant computed by pace2 + excitation energy (EE). (OK)
728	remmass = e_evapo;
729
730	G4double gamrem = std::sqrt(std::pow(pcorem,2)+std::pow(remmass,2))/remmass;
731	// assert(isnan(gamrem) == false);
732	G4double etrem = pcorem/remmass;
733
734	if(verboseLevel > 2) {
735	G4cout <<"5th Translab (no fission): Adding indices from " << 1 << " to " << volant->iv << G4endl;
736	}
737	nopart = varntp->ntrack - 1;
738	translab(gamrem,etrem,csrem,nopart,1);
739
740	// C End of the (FISSION - NO FISSION) condition (FF=1 or 0)
741	} //end if
742	// C ********************* FIN de l'EVAPO KHS ******************
743	}
744
745	// Evaporation code
746	void G4Abla::initEvapora()
747	{
748	// 37 C PROJECTILE AND TARGET PARAMETERS + CROSS SECTIONS
749	// 38 C COMMON /ABLAMAIN/ AP,ZP,AT,ZT,EAP,BETA,BMAXNUC,CRTOT,CRNUC,
750	// 39 C R_0,R_P,R_T, IMAX,IRNDM,PI,
751	// 40 C BFPRO,SNPRO,SPPRO,SHELL
752	// 41 C
753	// 42 C AP,ZP,AT,ZT - PROJECTILE AND TARGET MASSES
754	// 43 C EAP,BETA - BEAM ENERGY PER NUCLEON, V/C
755	// 44 C BMAXNUC - MAX. IMPACT PARAMETER FOR NUCL. REAC.
756	// 45 C CRTOT,CRNUC - TOTAL AND NUCLEAR REACTION CROSS SECTION
757	// 46 C R_0,R_P,R_T, - RADIUS PARAMETER, PROJECTILE+ TARGET RADII
758	// 47 C IMAX,IRNDM,PI - MAXIMUM NUMBER OF EVENTS, DUMMY, 3.141...
759	// 48 C BFPRO - FISSION BARRIER OF THE PROJECTILE
760	// 49 C SNPRO - NEUTRON SEPARATION ENERGY OF THE PROJECTILE
761	// 50 C SPPRO - PROTON " " " " "
762	// 51 C SHELL - GROUND STATE SHELL CORRECTION
763	// 52 C---------------------------------------------------------------------
764	// 53 C
765	// 54 C ENERGIES WIDTHS AND CROSS SECTIONS FOR EM EXCITATION
766	// 55 C COMMON /EMDPAR/ EGDR,EGQR,FWHMGDR,FWHMGQR,CREMDE1,CREMDE2,
767	// 56 C AE1,BE1,CE1,AE2,BE2,CE2,SR1,SR2,XR
768	// 57 C
769	// 58 C EGDR,EGQR - MEAN ENERGY OF GDR AND GQR
770	// 59 C FWHMGDR,FWHMGQR - FWHM OF GDR, GQR
771	// 60 C CREMDE1,CREMDE2 - EM CROSS SECTION FOR E1 AND E2
772	// 61 C AE1,BE1,CE1 - ARRAYS TO CALCULATE
773	// 62 C AE2,BE2,CE2 - THE EXCITATION ENERGY AFTER E.M. EXC.
774	// 63 C SR1,SR2,XR - WITH MONTE CARLO
775	// 64 C---------------------------------------------------------------------
776	// 65 C
777	// 66 C DEFORMATIONS AND G.S. SHELL EFFECTS
778	// 67 C COMMON /ECLD/ ECGNZ,ECFNZ,VGSLD,ALPHA
779	// 68 C
780	// 69 C ECGNZ - GROUND STATE SHELL CORR. FRLDM FOR A SPHERICAL G.S.
781	// 70 C ECFNZ - SHELL CORRECTION FOR THE SADDLE POINT (NOW: == 0)
782	// 71 C VGSLD - DIFFERENCE BETWEEN DEFORMED G.S. AND LDM VALUE
783	// 72 C ALPHA - ALPHA GROUND STATE DEFORMATION (THIS IS NOT BETA2!)
784	// 73 C BETA2 = SQRT(5/(4PI)) * ALPHA
785	// 74 C---------------------------------------------------------------------
786	// 75 C
787	// 76 C ARRAYS FOR EXCITATION ENERGY BY STATISTICAL HOLE ENERY MODEL
788	// 77 C COMMON /EENUC/ SHE, XHE
789	// 78 C
790	// 79 C SHE, XHE - ARRAYS TO CALCULATE THE EXC. ENERGY AFTER
791	// 80 C ABRASION BY THE STATISTICAL HOLE ENERGY MODEL
792	// 81 C---------------------------------------------------------------------
793	// 82 C
794	// 83 C G.S. SHELL EFFECT
795	// 84 C COMMON /EC2SUB/ ECNZ
796	// 85 C
797	// 86 C ECNZ G.S. SHELL EFFECT FOR THE MASSES (IDENTICAL TO ECGNZ)
798	// 87 C---------------------------------------------------------------------
799	// 88 C
800	// 89 C OPTIONS AND PARAMETERS FOR FISSION CHANNEL
801	// 90 C COMMON /FISS/ AKAP,BET,HOMEGA,KOEFF,IFIS,
802	// 91 C OPTSHP,OPTXFIS,OPTLES,OPTCOL
803	// 92 C
804	// 93 C AKAP - HBAR*2/(2 MN * R_0**2) = 10 MEV
805	// 94 C BET - REDUCED NUCLEAR FRICTION COEFFICIENT IN (1021 S-1)
806	// 95 C HOMEGA - CURVATURE OF THE FISSION BARRIER = 1 MEV
807	// 96 C KOEFF - COEFFICIENT FOR THE LD FISSION BARRIER == 1.0
808	// 97 C IFIS - 0/1 FISSION CHANNEL OFF/ON
809	// 98 C OPTSHP - INTEGER SWITCH FOR SHELL CORRECTION IN MASSES/ENERGY
810	// 99 C = 0 NO MICROSCOPIC CORRECTIONS IN MASSES AND ENERGY
811	// 100 C = 1 SHELL , NO PAIRING
812	// 101 C = 2 PAIRING, NO SHELL
813	// 102 C = 3 SHELL AND PAIRING
814	// 103 C OPTCOL - 0/1 COLLECTIVE ENHANCEMENT SWITCHED ON/OFF
815	// 104 C OPTXFIS- 0,1,2 FOR MYERS & SWIATECKI, DAHLINGER, ANDREYEV
816	// 105 C FISSILITY PARAMETER.
817	// 106 C OPTLES - CONSTANT TEMPERATURE LEVEL DENSITY FOR A,Z > TH-224
818	// 107 C OPTCOL - 0/1 COLLECTIVE ENHANCEMENT OFF/ON
819	// 108 C---------------------------------------------------------------------
820	// 109 C
821	// 110 C OPTIONS
822	// 111 C COMMON /OPT/ OPTEMD,OPTCHA,EEFAC
823	// 112 C
824	// 113 C OPTEMD - 0/1 NO EMD / INCL. EMD
825	// 114 C OPTCHA - 0/1 0 GDR / 1 HYPERGEOMETRICAL PREFRAGMENT-CHARGE-DIST.
826	// 115 C * RECOMMENDED IS OPTCHA = 1 *
827	// 116 C EEFAC - EXCITATION ENERGY FACTOR, 2.0 RECOMMENDED
828	// 117 C---------------------------------------------------------------------
829	// 118 C
830	// 119 C FISSION BARRIERS
831	// 120 C COMMON /FB/ EFA
832	// 121 C EFA - ARRAY OF FISSION BARRIERS
833	// 122 C---------------------------------------------------------------------
834	// 123 C
835	// 124 C p LEVEL DENSITY PARAMETERS
836	// 125 C COMMON /ALD/ AV,AS,AK,OPTAFAN
837	// 126 C AV,AS,AK - VOLUME,SURFACE,CURVATURE DEPENDENCE OF THE
838	// 127 C LEVEL DENSITY PARAMETER
839	// 128 C OPTAFAN - 0/1 AF/AN >=1 OR AF/AN ==1
840	// 129 C RECOMMENDED IS OPTAFAN = 0
841	// 130 C---------------------------------------------------------------------
842	// 131 C ____________________________________________________________________
843	// 132 C /
844	// 133 C / INITIALIZES PARAMETERS IN COMMON /ABRAMAIN/, /EMDPAR/, /ECLD/ ...
845	// 134 C / PROJECTILE PARAMETERS, EMD PARAMETERS, SHELL CORRECTION TABLES.
846	// 135 C / CALCULATES MAXIMUM IMPACT PARAMETER FOR NUCLEAR COLLISIONS AND
847	// 136 C / TOTAL GEOMETRICAL CROSS SECTION + EMD CROSS SECTIONS
848	// 137 C ____________________________________________________________________
849	// 138 C
850	// 139 C
851	// 201 C
852	// 202 C---------- SET INPUT VALUES
853	// 203 C
854	// 204 C *** INPUT FROM UNIT 10 IN THE FOLLOWING SEQUENCE !
855	// 205 C AP1 = INTEGER !
856	// 206 C ZP1 = INTEGER !
857	// 207 C AT1 = INTEGER !
858	// 208 C ZT1 = INTEGER !
859	// 209 C EAP = REAL !
860	// 210 C IMAX = INTEGER !
861	// 211 C IFIS = INTEGER SWITCH FOR FISSION
862	// 212 C OPTSHP = INTEGER SWITCH FOR SHELL CORRECTION IN MASSES/ENERGY
863	// 213 C =0 NO MICROSCOPIC CORRECTIONS IN MASSES AND ENERGY
864	// 214 C =1 SHELL , NO PAIRING CORRECTION
865	// 215 C =2 PAIRING, NO SHELL CORRECTION
866	// 216 C =3 SHELL AND PAIRING CORRECTION IN MASSES AND ENERGY
867	// 217 C OPTEMD =0,1 0 NO EMD, 1 INCL. EMD
868	// 218 C ELECTROMAGNETIC DISSOZIATION IS CALCULATED AS WELL.
869	// 219 C OPTCHA =0,1 0 GDR- , 1 HYPERGEOMETRICAL PREFRAGMENT-CHARGE-DIST.
870	// 220 C RECOMMENDED IS OPTCHA=1
871	// 221 C OPTCOL =0,1 COLLECTIVE ENHANCEMENT SWITCHED ON 1 OR OFF 0 IN DENSN
872	// 222 C OPTAFAN=0,1 SWITCH FOR AF/AN = 1 IN DENSNIV 0 AF/AN>1 1 AF/AN=1
873	// 223 C AKAP = REAL ALWAYS EQUALS 10
874	// 224 C BET = REAL REDUCED FRICTION COEFFICIENT / 10(+21) S(-1)
875	// 225 C HOMEGA = REAL CURVATURE / MEV RECOMMENDED = 1. MEV
876	// 226 C KOEFF = REAL COEFFICIENT FOR FISSION BARRIER
877	// 227 C OPTXFIS= INTEGER 0,1,2 FOR MYERS & SWIATECKI, DAHLINGER, ANDREYEV
878	// 228 C FISSILITY PARAMETER.
879	// 229 C EEFAC = REAL EMPIRICAL FACTOR FOR THE EXCITATION ENERGY
880	// 230 C RECOMMENDED 2.D0, STATISTICAL ABRASION MODELL 1.D0
881	// 231 C AV = REAL KOEFFICIENTS FOR CALCULATION OF A(TILDE)
882	// 232 C AS = REAL LEVEL DENSITY PARAMETER
883	// 233 C AK = REAL
884	// 234 C
885	// 235 C This following inputs will be initialized in the main through the
886	// 236 C common /ABLAMAIN/ (A.B.)
887	// 237
888
889	// switch-fission.1=on.0=off
890	fiss->ifis = 1;
891
892	// shell+pairing.0-1-2-3
893	fiss->optshp = 0;
894
895	// optemd =0,1 0 no emd, 1 incl. emd
896	opt->optemd = 1;
897	// read(10,*,iostat=io) dum(10),optcha
898	opt->optcha = 1;
899
900	// not.to.be.changed.(akap)
901	fiss->akap = 10.0;
902
903	// nuclear.viscosity.(beta)
904	fiss->bet = 1.5;
905
906	// potential-curvature
907	fiss->homega = 1.0;
908
909	// fission-barrier-coefficient
910	fiss->koeff = 1.;
911
912	//collective enhancement switched on 1 or off 0 in densn (qr=val or =1.)
913	fiss->optcol = 0;
914
915	// switch-for-low-energy-sys
916	fiss->optles = 0;
917
918	opt->eefac = 2.;
919
920	ald->optafan = 0;
921
922	ald->av = 0.073e0;
923	ald->as = 0.095e0;
924	ald->ak = 0.0e0;
925
926	if(verboseLevel > 3) {
927	G4cout <<"ifis " << fiss->ifis << G4endl;
928	G4cout <<"optshp " << fiss->optshp << G4endl;
929	G4cout <<"optemd " << opt->optemd << G4endl;
930	G4cout <<"optcha " << opt->optcha << G4endl;
931	G4cout <<"akap " << fiss->akap << G4endl;
932	G4cout <<"bet " << fiss->bet << G4endl;
933	G4cout <<"homega " << fiss->homega << G4endl;
934	G4cout <<"koeff " << fiss->koeff << G4endl;
935	G4cout <<"optcol " << fiss->optcol << G4endl;
936	G4cout <<"optles " << fiss->optles << G4endl;
937	G4cout <<"eefac " << opt->eefac << G4endl;
938	G4cout <<"optafan " << ald->optafan << G4endl;
939	G4cout <<"av " << ald->av << G4endl;
940	G4cout <<"as " << ald->as << G4endl;
941	G4cout <<"ak " << ald->ak << G4endl;
942	}
943	fiss->optxfis = 1;
944
945	G4InclAblaDataFile *dataInterface = new G4InclAblaDataFile();
946	if(dataInterface->readData() == true) {
947	if(verboseLevel > 0) {
948	G4cout <<"G4Abla: Datafiles read successfully." << G4endl;
949	}
950	}
951	else {
952	G4Exception("ERROR: Failed to read datafiles.");
953	}
954
955	for(int z = 0; z < 98; z++) { //do 30 z = 0,98,1
956	for(int n = 0; n < 154; n++) { //do 31 n = 0,153,1
957	ecld->ecfnz[n][z] = 0.e0;
958	ec2sub->ecnz[n][z] = dataInterface->getEcnz(n,z);
959	ecld->ecgnz[n][z] = dataInterface->getEcnz(n,z);
960	ecld->alpha[n][z] = dataInterface->getAlpha(n,z);
961	ecld->vgsld[n][z] = dataInterface->getVgsld(n,z);
962	}
963	}
964
965	for(int z = 0; z < 500; z++) {
966	for(int a = 0; a < 500; a++) {
967	pace->dm[z][a] = dataInterface->getPace2(z,a);
968	}
969	}
970
971	delete dataInterface;
972	}
973
974	void G4Abla::qrot(G4double z, G4double a, G4double bet, G4double sig, G4double u, G4double *qr)
975	{
976	G4double ucr = 10.0; // Critical energy for damping.
977	G4double dcr = 40.0; // Width of damping.
978	G4double ponq = 0.0, dn = 0.0, n = 0.0, dz = 0.0;
979
980	if(((std::fabs(bet)-1.15) < 0) \|\| ((std::fabs(bet)-1.15) == 0)) {
981	goto qrot10;
982	}
983
984	if((std::fabs(bet)-1.15) > 0) {
985	goto qrot11;
986	}
987
988	qrot10:
989	n = a - z;
990	dz = std::fabs(z - 82.0);
991	if (n > 104) {
992	dn = std::fabs(n-126.e0);
993	}
994	else {
995	dn = std::fabs(n - 82.0);
996	}
997
998	bet = 0.022 + 0.003dn + 0.005dz;
999
1000	sig = 25.0std::pow(bet,2) sig;
1001
1002	qrot11:
1003	ponq = (u - ucr)/dcr;
1004
1005	if (ponq > 700.0) {
1006	ponq = 700.0;
1007	}
1008	if (sig < 1.0) {
1009	sig = 1.0;
1010	}
1011	(qr) = 1.0/(1.0 + std::exp(ponq)) (sig - 1.0) + 1.0;
1012
1013	if ((*qr) < 1.0) {
1014	(*qr) = 1.0;
1015	}
1016
1017	return;
1018	}
1019
1020	void G4Abla::mglw(G4double a, G4double z, G4double *el)
1021	{
1022	// MODEL DE LA GOUTTE LIQUIDE DE C. F. WEIZSACKER.
1023	// USUALLY AN OBSOLETE OPTION
1024
1025	G4int a1 = 0, z1 = 0;
1026	G4double xv = 0.0, xs = 0.0, xc = 0.0, xa = 0.0;
1027
1028	a1 = idnint(a);
1029	z1 = idnint(z);
1030
1031	if ((a <= 0.01) \|\| (z < 0.01)) {
1032	(*el) = 1.0e38;
1033	}
1034	else {
1035	xv = -15.56*a;
1036	xs = 17.23*std::pow(a,(2.0/3.0));
1037
1038	if (a > 1.0) {
1039	xc = 0.7z(z-1.0)*std::pow((a-1.0),(-1.e0/3.e0));
1040	}
1041	else {
1042	xc = 0.0;
1043	}
1044	}
1045
1046	xa = 23.6(std::pow((a-2.0z),2)/a);
1047	(*el) = xv+xs+xc+xa;
1048	return;
1049	}
1050
1051	void G4Abla::mglms(G4double a, G4double z, G4int refopt4, G4double *el)
1052	{
1053	// USING FUNCTION EFLMAC(IA,IZ,0)
1054	//
1055	// REFOPT4 = 0 : WITHOUT MICROSCOPIC CORRECTIONS
1056	// REFOPT4 = 1 : WITH SHELL CORRECTION
1057	// REFOPT4 = 2 : WITH PAIRING CORRECTION
1058	// REFOPT4 = 3 : WITH SHELL- AND PAIRING CORRECTION
1059
1060	// 1839 C-----------------------------------------------------------------------
1061	// 1840 C A1 LOCAL MASS NUMBER (INTEGER VARIABLE OF A)
1062	// 1841 C Z1 LOCAL NUCLEAR CHARGE (INTEGER VARIABLE OF Z)
1063	// 1842 C REFOPT4 OPTION, SPECIFYING THE MASS FORMULA (SEE ABOVE)
1064	// 1843 C A MASS NUMBER
1065	// 1844 C Z NUCLEAR CHARGE
1066	// 1845 C DEL PAIRING CORRECTION
1067	// 1846 C EL BINDING ENERGY
1068	// 1847 C ECNZ( , ) TABLE OF SHELL CORRECTIONS
1069	// 1848 C-----------------------------------------------------------------------
1070	// 1849 C
1071	G4int a1 = idnint(a);
1072	G4int z1 = idnint(z);
1073
1074	if ( (a1 <= 0) \|\| (z1 <= 0) \|\| ((a1-z1) <= 0) ) { //then
1075	// modif pour rÃ©cupÃ©rer une masse p et n correcte:
1076	(*el) = 0.0;
1077	return;
1078	// goto mglms50;
1079	}
1080	else {
1081	// binding energy incl. pairing contr. is calculated from
1082	// function eflmac
1083	assert(a1 != 0);
1084	(*el) = eflmac(a1,z1,0,refopt4);
1085	// assert(isnan((*el)) == false);
1086	if (refopt4 > 0) {
1087	if (refopt4 != 2) {
1088	(el) = (el) + ec2sub->ecnz[a1-z1][z1];
1089	//(el) = (el) + ec2sub->ecnz[z1][a1-z1];
1090	//assert(isnan((*el)) == false);
1091	}
1092	}
1093	}
1094	return;
1095	}
1096
1097	G4double G4Abla::spdef(G4int a, G4int z, G4int optxfis)
1098	{
1099
1100	// INPUT: A,Z,OPTXFIS MASS AND CHARGE OF A NUCLEUS,
1101	// OPTION FOR FISSILITY
1102	// OUTPUT: SPDEF
1103
1104	// ALPHA2 SADDLE POINT DEF. COHEN&SWIATECKI ANN.PHYS. 22 (1963) 406
1105	// RANGING FROM FISSILITY X=0.30 TO X=1.00 IN STEPS OF 0.02
1106
1107	G4int index = 0;
1108	G4double x = 0.0, v = 0.0, dx = 0.0;
1109
1110	const G4int alpha2Size = 37;
1111	// The value 0.0 at alpha2[0] added by PK.
1112	G4double alpha2[alpha2Size] = {0.0, 2.5464e0, 2.4944e0, 2.4410e0, 2.3915e0, 2.3482e0,
1113	2.3014e0, 2.2479e0, 2.1982e0, 2.1432e0, 2.0807e0, 2.0142e0, 1.9419e0,
1114	1.8714e0, 1.8010e0, 1.7272e0, 1.6473e0, 1.5601e0, 1.4526e0, 1.3164e0,
1115	1.1391e0, 0.9662e0, 0.8295e0, 0.7231e0, 0.6360e0, 0.5615e0, 0.4953e0,
1116	0.4354e0, 0.3799e0, 0.3274e0, 0.2779e0, 0.2298e0, 0.1827e0, 0.1373e0,
1117	0.0901e0, 0.0430e0, 0.0000e0};
1118
1119	dx = 0.02;
1120	x = fissility(a,z,optxfis);
1121
1122	if (x > 1.0) {
1123	x = 1.0;
1124	}
1125
1126	if (x < 0.0) {
1127	x = 0.0;
1128	}
1129
1130	v = (x - 0.3)/dx + 1.0;
1131	index = idnint(v);
1132
1133	if (index < 1) {
1134	return(alpha2[1]); // alpha2[0] -> alpha2[1]
1135	}
1136
1137	if (index == 36) { //then // :::FIXME:: Possible off-by-one bug...
1138	return(alpha2[36]);
1139	}
1140	else {
1141	return(alpha2[index] + (alpha2[index+1] - alpha2[index]) / dx * ( x - (0.3e0 + dx*(index-1)))); //:::FIXME::: Possible off-by-one
1142	}
1143
1144	return alpha2[0]; // The algorithm is not supposed to reach this point.
1145	}
1146
1147	G4double G4Abla::fissility(int a,int z, int optxfis)
1148	{
1149	// CALCULATION OF FISSILITY PARAMETER
1150	//
1151	// INPUT: A,Z INTEGER MASS & CHARGE OF NUCLEUS
1152	// OPTXFIS = 0 : MYERS, SWIATECKI
1153	// 1 : DAHLINGER
1154	// 2 : ANDREYEV
1155
1156	G4double aa = 0.0, zz = 0.0, i = 0.0;
1157	G4double fissilityResult = 0.0;
1158
1159	aa = double(a);
1160	zz = double(z);
1161	i = double(a-2*z) / aa;
1162
1163	// myers & swiatecki droplet modell
1164	if (optxfis == 0) { //then
1165	fissilityResult = std::pow(zz,2) / aa /50.8830e0 / (1.0e0 - 1.7826e0 * std::pow(i,2));
1166	}
1167
1168	if (optxfis == 1) {
1169	// dahlinger fit:
1170	fissilityResult = std::pow(zz,2) / aa * std::pow((49.22e0(1.e0 - 0.3803e0std::pow(i,2) - 20.489e0*std::pow(i,4))),(-1));
1171	}
1172
1173	if (optxfis == 2) {
1174	// dubna fit:
1175	fissilityResult = std::pow(zz,2) / aa /(48.e0(1.e0 - 17.22e0std::pow(i,4)));
1176	}
1177
1178	return fissilityResult;
1179	}
1180
1181	void G4Abla::evapora(G4double zprf, G4double aprf, G4double ee, G4double jprf,
1182	G4double zf_par, G4double af_par, G4double *mtota_par,
1183	G4double pleva_par, G4double pxeva_par, G4double *pyeva_par,
1184	G4int ff_par, G4int inttype_par, G4int *inum_par)
1185	{
1186	G4double zf = (*zf_par);
1187	G4double af = (*af_par);
1188	G4double mtota = (*mtota_par);
1189	G4double pleva = (*pleva_par);
1190	G4double pxeva = (*pxeva_par);
1191	G4double pyeva = (*pyeva_par);
1192	G4int ff = (*ff_par);
1193	G4int inttype = (*inttype_par);
1194	G4int inum = (*inum_par);
1195
1196	// 533 C
1197	// 534 C INPUT:
1198	// 535 C
1199	// 536 C ZPRF, APRF, EE(EE IS MODIFIED!), JPRF
1200	// 537 C
1201	// 538 C PROJECTILE AND TARGET PARAMETERS + CROSS SECTIONS
1202	// 539 C COMMON /ABRAMAIN/ AP,ZP,AT,ZT,EAP,BETA,BMAXNUC,CRTOT,CRNUC,
1203	// 540 C R_0,R_P,R_T, IMAX,IRNDM,PI,
1204	// 541 C BFPRO,SNPRO,SPPRO,SHELL
1205	// 542 C
1206	// 543 C AP,ZP,AT,ZT - PROJECTILE AND TARGET MASSES
1207	// 544 C EAP,BETA - BEAM ENERGY PER NUCLEON, V/C
1208	// 545 C BMAXNUC - MAX. IMPACT PARAMETER FOR NUCL. REAC.
1209	// 546 C CRTOT,CRNUC - TOTAL AND NUCLEAR REACTION CROSS SECTION
1210	// 547 C R_0,R_P,R_T, - RADIUS PARAMETER, PROJECTILE+ TARGET RADII
1211	// 548 C IMAX,IRNDM,PI - MAXIMUM NUMBER OF EVENTS, DUMMY, 3.141...
1212	// 549 C BFPRO - FISSION BARRIER OF THE PROJECTILE
1213	// 550 C SNPRO - NEUTRON SEPARATION ENERGY OF THE PROJECTILE
1214	// 551 C SPPRO - PROTON " " " " "
1215	// 552 C SHELL - GROUND STATE SHELL CORRECTION
1216	// 553 C
1217	// 554 C---------------------------------------------------------------------
1218	// 555 C FISSION BARRIERS
1219	// 556 C COMMON /FB/ EFA
1220	// 557 C EFA - ARRAY OF FISSION BARRIERS
1221	// 558 C---------------------------------------------------------------------
1222	// 559 C OUTPUT:
1223	// 560 C ZF, AF, MTOTA, PLEVA, PTEVA, FF, INTTYPE, INUM
1224	// 561 C
1225	// 562 C ZF,AF - CHARGE AND MASS OF FINAL FRAGMENT AFTER EVAPORATION
1226	// 563 C MTOTA _ NUMBER OF EVAPORATED ALPHAS
1227	// 564 C PLEVA,PXEVA,PYEVA - MOMENTUM RECOIL BY EVAPORATION
1228	// 565 C INTTYPE - TYPE OF REACTION 0/1 NUCLEAR OR ELECTROMAGNETIC
1229	// 566 C FF - 0/1 NO FISSION / FISSION EVENT
1230	// 567 C INUM - EVENTNUMBER
1231	// 568 C ____________________________________________________________________
1232	// 569 C /
1233	// 570 C / CALCUL DE LA MASSE ET CHARGE FINALES D'UNE CHAINE D'EVAPORATION
1234	// 571 C /
1235	// 572 C / PROCEDURE FOR CALCULATING THE FINAL MASS AND CHARGE VALUES OF A
1236	// 573 C / SPECIFIC EVAPORATION CHAIN, STARTING POINT DEFINED BY (APRF, ZPRF,
1237	// 574 C / EE)
1238	// 575 C / On ajoute les 3 composantes de l'impulsion (PXEVA,PYEVA,PLEVA)
1239	// 576 C / (actuellement PTEVA n'est pas correct; mauvaise norme...)
1240	// 577 C /____________________________________________________________________
1241	// 578 C
1242	// 612 C
1243	// 613 C-----------------------------------------------------------------------
1244	// 614 C IRNDM DUMMY ARGUMENT FOR RANDOM-NUMBER FUNCTION
1245	// 615 C SORTIE LOCAL HELP VARIABLE TO END THE EVAPORATION CHAIN
1246	// 616 C ZF NUCLEAR CHARGE OF THE FRAGMENT
1247	// 617 C ZPRF NUCLEAR CHARGE OF THE PREFRAGMENT
1248	// 618 C AF MASS NUMBER OF THE FRAGMENT
1249	// 619 C APRF MASS NUMBER OF THE PREFRAGMENT
1250	// 620 C EPSILN ENERGY BURNED IN EACH EVAPORATION STEP
1251	// 621 C MALPHA LOCAL MASS CONTRIBUTION TO MTOTA IN EACH EVAPORATION
1252	// 622 C STEP
1253	// 623 C EE EXCITATION ENERGY (VARIABLE)
1254	// 624 C PROBP PROTON EMISSION PROBABILITY
1255	// 625 C PROBN NEUTRON EMISSION PROBABILITY
1256	// 626 C PROBA ALPHA-PARTICLE EMISSION PROBABILITY
1257	// 627 C PTOTL TOTAL EMISSION PROBABILITY
1258	// 628 C E LOWEST PARTICLE-THRESHOLD ENERGY
1259	// 629 C SN NEUTRON SEPARATION ENERGY
1260	// 630 C SBP PROTON SEPARATION ENERGY PLUS EFFECTIVE COULOMB
1261	// 631 C BARRIER
1262	// 632 C SBA ALPHA-PARTICLE SEPARATION ENERGY PLUS EFFECTIVE
1263	// 633 C COULOMB BARRIER
1264	// 634 C BP EFFECTIVE PROTON COULOMB BARRIER
1265	// 635 C BA EFFECTIVE ALPHA COULOMB BARRIER
1266	// 636 C MTOTA TOTAL MASS OF THE EVAPORATED ALPHA PARTICLES
1267	// 637 C X UNIFORM RANDOM NUMBER FOR NUCLEAR CHARGE
1268	// 638 C AMOINS LOCAL MASS NUMBER OF EVAPORATED PARTICLE
1269	// 639 C ZMOINS LOCAL NUCLEAR CHARGE OF EVAPORATED PARTICLE
1270	// 640 C ECP KINETIC ENERGY OF PROTON WITHOUT COULOMB
1271	// 641 C REPULSION
1272	// 642 C ECN KINETIC ENERGY OF NEUTRON
1273	// 643 C ECA KINETIC ENERGY OF ALPHA PARTICLE WITHOUT COULOMB
1274	// 644 C REPULSION
1275	// 645 C PLEVA TRANSVERSAL RECOIL MOMENTUM OF EVAPORATION
1276	// 646 C PTEVA LONGITUDINAL RECOIL MOMENTUM OF EVAPORATION
1277	// 647 C FF FISSION FLAG
1278	// 648 C INTTYPE INTERACTION TYPE FLAG
1279	// 649 C RNDX RECOIL MOMENTUM IN X-DIRECTION IN A SINGLE STEP
1280	// 650 C RNDY RECOIL MOMENTUM IN Y-DIRECTION IN A SINGLE STEP
1281	// 651 C RNDZ RECOIL MOMENTUM IN Z-DIRECTION IN A SINGLE STEP
1282	// 652 C RNDN NORMALIZATION OF RECOIL MOMENTUM FOR EACH STEP
1283	// 653 C-----------------------------------------------------------------------
1284	// 654 C
1285	// 655 SAVE
1286	// SAVE -> static
1287
1288	static G4int sortie = 0;
1289	static G4double epsiln = 0.0, probp = 0.0, probn = 0.0, proba = 0.0, ptotl = 0.0, e = 0.0;
1290	static G4double sn = 0.0, sbp = 0.0, sba = 0.0, x = 0.0, amoins = 0.0, zmoins = 0.0;
1291	G4double ecn = 0.0, ecp = 0.0,eca = 0.0, bp = 0.0, ba = 0.0;
1292	static G4double pteva = 0.0;
1293
1294	static G4int itest = 0;
1295	static G4double probf = 0.0;
1296
1297	static G4int k = 0, j = 0, il = 0;
1298
1299	static G4double ctet1 = 0.0, stet1 = 0.0, phi1 = 0.0;
1300	static G4double sbfis = 0.0, rnd = 0.0;
1301	static G4double selmax = 0.0;
1302	static G4double segs = 0.0;
1303	static G4double ef = 0.0;
1304	static G4int irndm = 0;
1305
1306	static G4double pc = 0.0, malpha = 0.0;
1307
1308	zf = zprf;
1309	af = aprf;
1310	pleva = 0.0;
1311	pteva = 0.0;
1312	pxeva = 0.0;
1313	pyeva = 0.0;
1314
1315	sortie = 0;
1316	ff = 0;
1317
1318	itest = 0;
1319	if (itest == 1) {
1320	G4cout << "***************************" << G4endl;
1321	}
1322
1323	evapora10:
1324
1325	if (itest == 1) {
1326	G4cout <<"------zf,af,ee------" << idnint(zf) << "," << idnint(af) << "," << ee << G4endl;
1327	}
1328
1329	// calculation of the probabilities for the different decay channels
1330	// plus separation energies and kinetic energies of the particles
1331	direct(zf,af,ee,jprf,&probp,&probn,&proba,&probf,&ptotl,
1332	&sn,&sbp,&sba,&ecn,&ecp,&eca,&bp,&ba,inttype,inum,itest); //:::FIXME::: Call
1333	// assert(isnan(proba) == false);
1334	// assert(isnan(probp) == false);
1335	// assert(isnan(probn) == false);
1336	// assert(isnan(probf) == false);
1337	assert((eca+ba) >= 0);
1338	assert((ecp+bp) >= 0);
1339	// assert(isnan(ecp) == false);
1340	// assert(isnan(ecn) == false);
1341	// assert(isnan(bp) == false);
1342	// assert(isnan(ba) == false);
1343	k = idnint(zf);
1344	j = idnint(af-zf);
1345
1346	// now ef is calculated from efa that depends on the subroutine
1347	// barfit which takes into account the modification on the ang. mom.
1348	// jb mvr 6-aug-1999
1349	// note *** shell correction! (ecgnz) jb mvr 20-7-1999
1350	il = idnint(jprf);
1351	barfit(k,k+j,il,&sbfis,&segs,&selmax);
1352	// assert(isnan(sbfis) == false);
1353
1354	if ((fiss->optshp == 1) \|\| (fiss->optshp == 3)) { //then
1355	// fb->efa[k][j] = G4double(sbfis) - ecld->ecgnz[j][k];
1356	fb->efa[j][k] = G4double(sbfis) - ecld->ecgnz[j][k];
1357	}
1358	else {
1359	fb->efa[j][k] = G4double(sbfis);
1360	// fb->efa[j][k] = G4double(sbfis);
1361	} //end if
1362	ef = fb->efa[j][k];
1363	// ef = fb->efa[j][k];
1364	// assert(isnan(fb->efa[j][k]) == false);
1365	// here the final steps of the evaporation are calculated
1366	if ((sortie == 1) \|\| (ptotl == 0.e0)) {
1367	e = dmin1(sn,sbp,sba);
1368	if (e > 1.0e30) {
1369	if(verboseLevel > 2) {
1370	G4cout <<"erreur a la sortie evapora,e>1.e30,af=" << af <<" zf=" << zf << G4endl;
1371	}
1372	}
1373	if (zf <= 6.0) {
1374	goto evapora100;
1375	}
1376	if (e < 0.0) {
1377	if (sn == e) {
1378	af = af - 1.e0;
1379	}
1380	else if (sbp == e) {
1381	af = af - 1.0;
1382	zf = zf - 1.0;
1383	}
1384	else if (sba == e) {
1385	af = af - 4.0;
1386	zf = zf - 2.0;
1387	}
1388	if (af < 2.5) {
1389	goto evapora100;
1390	}
1391	goto evapora10;
1392	}
1393	goto evapora100;
1394	}
1395	irndm = irndm + 1;
1396
1397	// here the normal evaporation cascade starts
1398
1399	// random number for the evaporation
1400	// x = double(Rndm(irndm))*ptotl;
1401	x = double(haz(1))*ptotl;
1402
1403	// G4cout <<"proba = " << proba << G4endl;
1404	// G4cout <<"probp = " << probp << G4endl;
1405	// G4cout <<"probn = " << probn << G4endl;
1406	// G4cout <<"probf = " << probf << G4endl;
1407
1408	itest = 0;
1409	if (x < proba) {
1410	// alpha evaporation
1411	if (itest == 1) {
1412	G4cout <<"< alpha evaporation >" << G4endl;
1413	}
1414	amoins = 4.0;
1415	zmoins = 2.0;
1416	epsiln = sba + eca;
1417	assert((std::pow((1.0 + (eca+ba)/3.72834e3),2) - 1.0) >= 0);
1418	pc = std::sqrt(std::pow((1.0 + (eca+ba)/3.72834e3),2) - 1.0) * 3.72834e3;
1419	// assert(isnan(pc) == false);
1420	malpha = 4.0;
1421
1422	// volant:
1423	volant->iv = volant->iv + 1;
1424	volant->acv[volant->iv] = 4.;
1425	volant->zpcv[volant->iv] = 2.;
1426	volant->pcv[volant->iv] = pc;
1427	}
1428	else if (x < proba+probp) {
1429	// proton evaporation
1430	if (itest == 1) {
1431	G4cout <<"< proton evaporation >" << G4endl;
1432	}
1433	amoins = 1.0;
1434	zmoins = 1.0;
1435	epsiln = sbp + ecp;
1436	assert((std::pow((1.0 + (ecp + bp)/9.3827e2),2) - 1.0) >= 0);
1437	pc = std::sqrt(std::pow((1.0 + (ecp + bp)/9.3827e2),2) - 1.0) * 9.3827e2;
1438	// assert(isnan(pc) == false);
1439	malpha = 0.0;
1440	// volant:
1441	volant->iv = volant->iv + 1;
1442	volant->acv[volant->iv] = 1.0;
1443	volant->zpcv[volant->iv] = 1.;
1444	volant->pcv[volant->iv] = pc;
1445	}
1446	else if (x < proba+probp+probn) {
1447	// neutron evaporation
1448	if (itest == 1) {
1449	G4cout <<"< neutron evaporation >" << G4endl;
1450	}
1451	amoins = 1.0;
1452	zmoins = 0.0;
1453	epsiln = sn + ecn;
1454	assert((std::pow((1.0 + (ecn)/9.3956e2),2) - 1.0) >= 0);
1455	pc = std::sqrt(std::pow((1.0 + (ecn)/9.3956e2),2) - 1.0) * 9.3956e2;
1456	// assert(isnan(pc) == false);
1457	malpha = 0.0;
1458
1459	// volant:
1460	volant->iv = volant->iv + 1;
1461	volant->acv[volant->iv] = 1.;
1462	volant->zpcv[volant->iv] = 0.;
1463	volant->pcv[volant->iv] = pc;
1464	}
1465	else {
1466	// fission
1467	// in case of fission-events the fragment nucleus is the mother nucleus
1468	// before fission occurs with excitation energy above the fis.- barrier.
1469	// fission fragment mass distribution is calulated in subroutine fisdis
1470	if (itest == 1) {
1471	G4cout <<"< fission >" << G4endl;
1472	}
1473	amoins = 0.0;
1474	zmoins = 0.0;
1475	epsiln = ef;
1476
1477	malpha = 0.0;
1478	pc = 0.0;
1479	ff = 1;
1480	// ff = 0; // For testing, allows to disable fission!
1481	}
1482
1483	if (itest == 1) {
1484	G4cout <<"sn,sbp,sba,ef" << sn << "," << sbp << "," << sba <<"," << ef << G4endl;
1485	G4cout <<"probn,probp,proba,probf,ptotl " <<","<< probn <<","<< probp <<","<< proba <<","<< probf <<","<< ptotl << G4endl;
1486	}
1487
1488	// calculation of the daughter nucleus
1489	af = af - amoins;
1490	zf = zf - zmoins;
1491	ee = ee - epsiln;
1492	if (ee <= 0.01) {
1493	ee = 0.01;
1494	}
1495	mtota = mtota + malpha;
1496
1497	if(ff == 0) {
1498	standardRandom(&rnd,&(hazard->igraine[8]));
1499	ctet1 = 2.0*rnd - 1.0;
1500	standardRandom(&rnd,&(hazard->igraine[4]));
1501	phi1 = rnd2.03.141592654;
1502	stet1 = std::sqrt(1.0 - std::pow(ctet1,2));
1503	// assert(isnan(stet1) == false);
1504	volant->xcv[volant->iv] = stet1*std::cos(phi1);
1505	volant->ycv[volant->iv] = stet1*std::sin(phi1);
1506	volant->zcv[volant->iv] = ctet1;
1507	pxeva = pxeva - pc * volant->xcv[volant->iv];
1508	pyeva = pyeva - pc * volant->ycv[volant->iv];
1509	pleva = pleva - pc * ctet1;
1510	// assert(isnan(pleva) == false);
1511	}
1512
1513	// condition for end of evaporation
1514	if ((af < 2.5) \|\| (ff == 1)) {
1515	goto evapora100;
1516	}
1517	goto evapora10;
1518
1519	evapora100:
1520	(*zf_par) = zf;
1521	(*af_par) = af;
1522	(*mtota_par) = mtota;
1523	(*pleva_par) = pleva;
1524	(*pxeva_par) = pxeva;
1525	(*pyeva_par) = pyeva;
1526	(*ff_par) = ff;
1527	(*inttype_par) = inttype;
1528	(*inum_par) = inum;
1529
1530	return;
1531	}
1532
1533	void G4Abla::direct(G4double zprf, G4double a, G4double ee, G4double jprf,
1534	G4double probp_par, G4double probn_par, G4double *proba_par,
1535	G4double probf_par, G4double ptotl_par, G4double *sn_par,
1536	G4double sbp_par, G4double sba_par, G4double *ecn_par,
1537	G4double ecp_par,G4double eca_par, G4double *bp_par,
1538	G4double *ba_par, G4int inttype, G4int inum, G4int itest)
1539	{
1540	G4int dummy0 = 0;
1541
1542	G4double probp = (*probp_par);
1543	G4double probn = (*probn_par);
1544	G4double proba = (*proba_par);
1545	G4double probf = (*probf_par);
1546	G4double ptotl = (*ptotl_par);
1547	G4double sn = (*sn_par);
1548	G4double sbp = (*sbp_par);
1549	G4double sba = (*sba_par);
1550	G4double ecn = (*ecn_par);
1551	G4double ecp = (*ecp_par);
1552	G4double eca = (*eca_par);
1553	G4double bp = (*bp_par);
1554	G4double ba = (*ba_par);
1555
1556	// CALCULATION OF PARTICLE-EMISSION PROBABILITIES & FISSION /
1557	// BASED ON THE SIMPLIFIED FORMULAS FOR THE DECAY WIDTH BY /
1558	// MORETTO, ROCHESTER MEETING TO AVOID COMPUTING TIME /
1559	// INTENSIVE INTEGRATION OF THE LEVEL DENSITIES /
1560	// USES EFFECTIVE COULOMB BARRIERS AND AN AVERAGE KINETIC ENERGY/
1561	// OF THE EVAPORATED PARTICLES /
1562	// COLLECTIVE ENHANCMENT OF THE LEVEL DENSITY IS INCLUDED /
1563	// DYNAMICAL HINDRANCE OF FISSION IS INCLUDED BY A STEP FUNCTION/
1564	// APPROXIMATION. SEE A.R. JUNGHANS DIPLOMA THESIS /
1565	// SHELL AND PAIRING STRUCTURES IN THE LEVEL DENSITY IS INCLUDED/
1566
1567	// INPUT:
1568	// ZPRF,A,EE CHARGE, MASS, EXCITATION ENERGY OF COMPOUND
1569	// NUCLEUS
1570	// JPRF ROOT-MEAN-SQUARED ANGULAR MOMENTUM
1571
1572	// DEFORMATIONS AND G.S. SHELL EFFECTS
1573	// COMMON /ECLD/ ECGNZ,ECFNZ,VGSLD,ALPHA
1574
1575	// ECGNZ - GROUND STATE SHELL CORR. FRLDM FOR A SPHERICAL G.S.
1576	// ECFNZ - SHELL CORRECTION FOR THE SADDLE POINT (NOW: == 0)
1577	// VGSLD - DIFFERENCE BETWEEN DEFORMED G.S. AND LDM VALUE
1578	// ALPHA - ALPHA GROUND STATE DEFORMATION (THIS IS NOT BETA2!)
1579	// BETA2 = SQRT((4PI)/5) * ALPHA
1580
1581	//OPTIONS AND PARAMETERS FOR FISSION CHANNEL
1582	//COMMON /FISS/ AKAP,BET,HOMEGA,KOEFF,IFIS,
1583	// OPTSHP,OPTXFIS,OPTLES,OPTCOL
1584	//
1585	// AKAP - HBAR*2/(2 MN * R_0**2) = 10 MEV, R_0 = 1.4 FM
1586	// BET - REDUCED NUCLEAR FRICTION COEFFICIENT IN (1021 S-1)
1587	// HOMEGA - CURVATURE OF THE FISSION BARRIER = 1 MEV
1588	// KOEFF - COEFFICIENT FOR THE LD FISSION BARRIER == 1.0
1589	// IFIS - 0/1 FISSION CHANNEL OFF/ON
1590	// OPTSHP - INTEGER SWITCH FOR SHELL CORRECTION IN MASSES/ENERGY
1591	// = 0 NO MICROSCOPIC CORRECTIONS IN MASSES AND ENERGY
1592	// = 1 SHELL , NO PAIRING
1593	// = 2 PAIRING, NO SHELL
1594	// = 3 SHELL AND PAIRING
1595	// OPTCOL - 0/1 COLLECTIVE ENHANCEMENT SWITCHED ON/OFF
1596	// OPTXFIS- 0,1,2 FOR MYERS & SWIATECKI, DAHLINGER, ANDREYEV
1597	// FISSILITY PARAMETER.
1598	// OPTLES - CONSTANT TEMPERATURE LEVEL DENSITY FOR A,Z > TH-224
1599	// OPTCOL - 0/1 COLLECTIVE ENHANCEMENT OFF/ON
1600
1601	// LEVEL DENSITY PARAMETERS
1602	// COMMON /ALD/ AV,AS,AK,OPTAFAN
1603	// AV,AS,AK - VOLUME,SURFACE,CURVATURE DEPENDENCE OF THE
1604	// LEVEL DENSITY PARAMETER
1605	// OPTAFAN - 0/1 AF/AN >=1 OR AF/AN ==1
1606	// RECOMMENDED IS OPTAFAN = 0
1607
1608	// FISSION BARRIERS
1609	// COMMON /FB/ EFA
1610	// EFA - ARRAY OF FISSION BARRIERS
1611
1612
1613	// OUTPUT: PROBN,PROBP,PROBA,PROBF,PTOTL:
1614	// - EMISSION PROBABILITIES FOR N EUTRON, P ROTON, A LPHA
1615	// PARTICLES, F ISSION AND NORMALISATION
1616	// SN,SBP,SBA: SEPARATION ENERGIES N P A
1617	// INCLUDING EFFECTIVE BARRIERS
1618	// ECN,ECP,ECA,BP,BA
1619	// - AVERAGE KINETIC ENERGIES (2*T) AND EFFECTIVE BARRIERS
1620
1621	static G4double bk = 0.0;
1622	static G4int afp = 0;
1623	static G4double at = 0.0;
1624	static G4double bs = 0.0;
1625	static G4double bshell = 0.0;
1626	static G4double cf = 0.0;
1627	static G4double dconst = 0.0;
1628	static G4double defbet = 0.0;
1629	static G4double denomi = 0.0;
1630	static G4double densa = 0.0;
1631	static G4double densf = 0.0;
1632	static G4double densg = 0.0;
1633	static G4double densn = 0.0;
1634	static G4double densp = 0.0;
1635	static G4double edyn = 0.0;
1636	static G4double eer = 0.0;
1637	static G4double ef = 0.0;
1638	static G4double ft = 0.0;
1639	static G4double ga = 0.0;
1640	static G4double gf = 0.0;
1641	static G4double gn = 0.0;
1642	static G4double gngf = 0.0;
1643	static G4double gp = 0.0;
1644	static G4double gsum = 0.0;
1645	static G4double hbar = 6.582122e-22; // = 0.0;
1646	static G4double iflag = 0.0;
1647	static G4int il = 0;
1648	static G4int imaxwell = 0;
1649	static G4int in = 0;
1650	static G4int iz = 0;
1651	static G4int j = 0;
1652	static G4int k = 0;
1653	static G4double ma1z = 0.0;
1654	static G4double ma1z1 = 0.0;
1655	static G4double ma4z2 = 0.0;
1656	static G4double maz = 0.0;
1657	static G4double nprf = 0.0;
1658	static G4double nt = 0.0;
1659	static G4double parc = 0.0;
1660	static G4double pi = 3.14159265;
1661	static G4double pt = 0.0;
1662	static G4double ra = 0.0;
1663	static G4double rat = 0.0;
1664	static G4double refmod = 0.0;
1665	static G4double rf = 0.0;
1666	static G4double rn = 0.0;
1667	static G4double rnd = 0.0;
1668	static G4double rnt = 0.0;
1669	static G4double rp = 0.0;
1670	static G4double rpt = 0.0;
1671	static G4double sa = 0.0;
1672	static G4double sbf = 0.0;
1673	static G4double sbfis = 0.0;
1674	static G4double segs = 0.0;
1675	static G4double selmax = 0.0;
1676	static G4double sp = 0.0;
1677	static G4double tauc = 0.0;
1678	static G4double tconst = 0.0;
1679	static G4double temp = 0.0;
1680	static G4double ts1 = 0.0;
1681	static G4double tsum = 0.0;
1682	static G4double wf = 0.0;
1683	static G4double wfex = 0.0;
1684	static G4double xx = 0.0;
1685	static G4double y = 0.0;
1686
1687	imaxwell = 1;
1688	inttype = 0;
1689
1690	// limiting of excitation energy where fission occurs
1691	// Note, this is not the dynamical hindrance (see end of routine)
1692	edyn = 1000.0;
1693
1694	// no limit if statistical model is calculated.
1695	if (fiss->bet <= 1.0e-16) {
1696	edyn = 10000.0;
1697	}
1698
1699	// just a change of name until the end of this subroutine
1700	eer = ee;
1701	if (inum == 1) {
1702	ilast = 1;
1703	}
1704
1705	// calculation of masses
1706	// refmod = 1 ==> myers,swiatecki model
1707	// refmod = 0 ==> weizsaecker model
1708	refmod = 1; // Default = 1
1709
1710	if (refmod == 1) {
1711	mglms(a,zprf,fiss->optshp,&maz);
1712	mglms(a-1.0,zprf,fiss->optshp,&ma1z);
1713	mglms(a-1.0,zprf-1.0,fiss->optshp,&ma1z1);
1714	mglms(a-4.0,zprf-2.0,fiss->optshp,&ma4z2);
1715	}
1716	else {
1717	mglw(a,zprf,&maz);
1718	mglw(a-1.0,zprf,&ma1z);
1719	mglw(a-1.0,zprf-1.0,&ma1z1);
1720	mglw(a-4.0,zprf-2.0,&ma4z2);
1721	}
1722	// assert(isnan(maz) == false);
1723	// assert(isnan(ma1z) == false);
1724	// assert(isnan(ma1z1) == false);
1725	// assert(isnan(ma4z2) == false);
1726
1727	// separation energies and effective barriers
1728	sn = ma1z - maz;
1729	sp = ma1z1 - maz;
1730	sa = ma4z2 - maz - 28.29688;
1731	if (zprf < 1.0e0) {
1732	sbp = 1.0e75;
1733	goto direct30;
1734	}
1735
1736	// parameterisation gaimard:
1737	// bp = 1.44(zprf-1.d0)/(1.22std::pow((a - 1.0),(1.0/3.0))+5.6)
1738	// parameterisation khs (12-99)
1739	bp = 1.44(zprf - 1.0)/(2.1std::pow((a - 1.0),(1.0/3.0)) + 0.0);
1740
1741	sbp = sp + bp;
1742	// assert(isnan(sbp) == false);
1743	// assert(isinf(sbp) == false);
1744	if (a-4.0 <= 0.0) {
1745	sba = 1.0e+75;
1746	goto direct30;
1747	}
1748
1749	// new effective barrier for alpha evaporation d=6.1: khs
1750	// ba = 2.88d0(zprf-2.d0)/(1.22d0(a-4.d0)**(1.d0/3.d0)+6.1d0)
1751	// parametrisation khs (12-99)
1752	ba = 2.88(zprf - 2.0)/(2.2std::pow((a - 4.0),(1.0/3.0)) + 0.0);
1753
1754	sba = sa + ba;
1755	// assert(isnan(sba) == false);
1756	// assert(isinf(sba) == false);
1757	direct30:
1758
1759	// calculation of surface and curvature integrals needed to
1760	// to calculate the level density parameter (in densniv)
1761	if (fiss->ifis > 0) {
1762	k = idnint(zprf);
1763	j = idnint(a - zprf);
1764
1765	// now ef is calculated from efa that depends on the subroutine
1766	// barfit which takes into account the modification on the ang. mom.
1767	// jb mvr 6-aug-1999
1768	// note *** shell correction! (ecgnz) jb mvr 20-7-1999
1769	il = idnint(jprf);
1770	barfit(k,k+j,il,&sbfis,&segs,&selmax);
1771	// assert(isnan(sbfis) == false);
1772	if ((fiss->optshp == 1) \|\| (fiss->optshp == 3)) {
1773	// fb->efa[k][j] = G4double(sbfis) - ecld->ecgnz[j][k];
1774	// fb->efa[j][k] = G4double(sbfis) - ecld->ecgnz[j][k];
1775	fb->efa[j][k] = double(sbfis) - ecld->ecgnz[j][k];
1776	}
1777	else {
1778	// fb->efa[k][j] = G4double(sbfis);
1779	fb->efa[j][k] = double(sbfis);
1780	}
1781	// ef = fb->efa[k][j];
1782	ef = fb->efa[j][k];
1783
1784	// to avoid negative values for impossible nuclei
1785	// the fission barrier is set to zero if smaller than zero.
1786	// if (fb->efa[k][j] < 0.0) {
1787	// fb->efa[k][j] = 0.0;
1788	// }
1789	if (fb->efa[j][k] < 0.0) {
1790	if(verboseLevel > 2) {
1791	G4cout <<"Setting fission barrier to 0" << G4endl;
1792	}
1793	fb->efa[j][k] = 0.0;
1794	}
1795	// assert(isnan(fb->efa[j][k]) == false);
1796
1797	// factor with jprf should be 0.0025d0 - 0.01d0 for
1798	// approximate influence of ang. momentum on bfis a.j. 22.07.96
1799	// 0.0 means no angular momentum
1800
1801	if (ef < 0.0) {
1802	ef = 0.0;
1803	}
1804	xx = fissility((k+j),k,fiss->optxfis);
1805	// assert(isnan(xx) == false);
1806	// assert(isinf(xx) == false);
1807
1808	y = 1.00 - xx;
1809	if (y < 0.0) {
1810	y = 0.0;
1811	}
1812	if (y > 1.0) {
1813	y = 1.0;
1814	}
1815	bs = bipol(1,y);
1816	// assert(isnan(bs) == false);
1817	// assert(isinf(bs) == false);
1818	bk = bipol(2,y);
1819	// assert(isnan(bk) == false);
1820	// assert(isinf(bk) == false);
1821	}
1822	else {
1823	ef = 1.0e40;
1824	bs = 1.0;
1825	bk = 1.0;
1826	}
1827	sbf = ee - ef;
1828
1829	afp = idnint(a);
1830	iz = idnint(zprf);
1831	in = afp - iz;
1832	bshell = ecld->ecfnz[in][iz];
1833	// assert(isnan(bshell) == false);
1834
1835	// ld saddle point deformation
1836	// here: beta2 = std::sqrt(5/(4pi)) * alpha2
1837
1838	// for the ground state def. 1.5d0 should be used
1839	// because this was just the factor to produce the
1840	// alpha-deformation table 1.5d0 should be used
1841	// a.r.j. 6.8.97
1842	defbet = 1.58533e0 * spdef(idnint(a),idnint(zprf),fiss->optxfis);
1843	// assert(isnan(defbet) == false);
1844
1845	// level density and temperature at the saddle point
1846	// G4cout <<"a = " << a << G4endl;
1847	// G4cout <<"zprf = " << zprf << G4endl;
1848	// G4cout <<"ee = " << ee << G4endl;
1849	// G4cout <<"ef = " << ef << G4endl;
1850	// G4cout <<"bshell = " << bshell << G4endl;
1851	// G4cout <<"bs = " << bs << G4endl;
1852	// G4cout <<"bk = " << bk << G4endl;
1853	// G4cout <<"defbet = " << defbet << G4endl;
1854	densniv(a,zprf,ee,ef,&densf,bshell,bs,bk,&temp,int(fiss->optshp),int(fiss->optcol),defbet);
1855	// G4cout <<"densf = " << densf << G4endl;
1856	// G4cout <<"temp = " << temp << G4endl;
1857	// assert(isnan(densf) == false);
1858	// assert(isnan(temp) == false);
1859	// assert(temp != 0);
1860	ft = temp;
1861	if (iz >= 2) {
1862	bshell = ecld->ecgnz[in][iz-1] - ecld->vgsld[in][iz-1];
1863	defbet = 1.5 * (ecld->alpha[in][iz-1]);
1864
1865	// level density and temperature in the proton daughter
1866	densniv(a-1.0,zprf-1.0e0,ee,sbp,&densp, bshell,1.e0,1.e0,&temp,int(fiss->optshp),int(fiss->optcol),defbet);
1867	assert(temp >= 0);
1868	// assert(isnan(temp) == false);
1869	pt = temp;
1870	if (imaxwell == 1) {
1871	// valentina - random kinetic energy in a maxwelliam distribution
1872	// modif juin/2002 a.b. c.v. for light targets; limit on the energy
1873	// from the maxwell distribution.
1874	rpt = pt;
1875	ecp = 2.0 * pt;
1876	if(rpt <= 1.0e-3) {
1877	goto direct2914;
1878	}
1879	iflag = 0;
1880	direct1914:
1881	ecp = fmaxhaz(rpt);
1882	iflag = iflag + 1;
1883	if(iflag >= 10) {
1884	standardRandom(&rnd,&(hazard->igraine[5]));
1885	ecp=std::sqrt(rnd)*(eer-sbp);
1886	// assert(isnan(ecp) == false);
1887	goto direct2914;
1888	}
1889	if((ecp+sbp) > eer) {
1890	goto direct1914;
1891	}
1892	}
1893	else {
1894	ecp = 2.0 * pt;
1895	}
1896
1897	direct2914:
1898	dummy0 = 0;
1899	// G4cout <<""<<G4endl;
1900	}
1901	else {
1902	densp = 0.0;
1903	ecp = 0.0;
1904	pt = 0.0;
1905	}
1906
1907	if (in >= 2) {
1908	bshell = ecld->ecgnz[in-1][iz] - ecld->vgsld[in-1][iz];
1909	defbet = 1.5e0 * (ecld->alpha[in-1][iz]);
1910
1911	// level density and temperature in the neutron daughter
1912	densniv(a-1.0,zprf,ee,sn,&densn,bshell, 1.e0,1.e0,&temp,int(fiss->optshp),int(fiss->optcol),defbet);
1913	nt = temp;
1914
1915	if (imaxwell == 1) {
1916	// valentina - random kinetic energy in a maxwelliam distribution
1917	// modif juin/2002 a.b. c.v. for light targets; limit on the energy
1918	// from the maxwell distribution.
1919	rnt = nt;
1920	ecn = 2.0 * nt;
1921	if(rnt <= 1.e-3) {
1922	goto direct2915;
1923	}
1924
1925	iflag=0;
1926
1927	ecn = fmaxhaz(rnt);
1928	iflag=iflag+1;
1929	if(iflag >= 10) {
1930	standardRandom(&rnd,&(hazard->igraine[6]));
1931	ecn = std::sqrt(rnd)*(eer-sn);
1932	// assert(isnan(ecn) == false);
1933	goto direct2915;
1934	}
1935	// if((ecn+sn) > eer) {
1936	// goto direct1915;
1937	// }
1938	// else {
1939	// ecn = 2.e0 * nt;
1940	// }
1941	if((ecn + sn) <= eer) {
1942	ecn = 2.0 * nt;
1943	}
1944	direct2915:
1945	dummy0 = 0;
1946	// G4cout <<"" <<G4endl;
1947	}
1948	}
1949	else {
1950	densn = 0.0;
1951	ecn = 0.0;
1952	nt = 0.0;
1953	}
1954
1955	if ((in >= 3) && (iz >= 3)) {
1956	bshell = ecld->ecgnz[in-2][iz-2] - ecld->vgsld[in-2][iz-2];
1957	defbet = 1.5 * (ecld->alpha[in-2][iz-2]);
1958
1959	// level density and temperature in the alpha daughter
1960	densniv(a-4.0,zprf-2.0e0,ee,sba,&densa,bshell,1.e0,1.e0,&temp,int(fiss->optshp),int(fiss->optcol),defbet);
1961
1962	// valentina - random kinetic energy in a maxwelliam distribution
1963	at = temp;
1964	if (imaxwell == 1) {
1965	// modif juin/2002 a.b. c.v. for light targets; limit on the energy
1966	// from the maxwell distribution.
1967	rat = at;
1968	eca= 2.e0 * at;
1969	if(rat <= 1.e-3) {
1970	goto direct2916;
1971	}
1972	iflag=0;
1973	direct1916:
1974	eca = fmaxhaz(rat);
1975	iflag=iflag+1;
1976	if(iflag >= 10) {
1977	standardRandom(&rnd,&(hazard->igraine[7]));
1978	eca=std::sqrt(rnd)*(eer-sba);
1979	// assert(isnan(eca) == false);
1980	goto direct2916;
1981	}
1982	if((eca+sba) > eer) {
1983	goto direct1916;
1984	}
1985	else {
1986	eca = 2.0 * at;
1987	}
1988	direct2916:
1989	dummy0 = 0;
1990	// G4cout <<"" << G4endl;
1991	}
1992	else {
1993	densa = 0.0;
1994	eca = 0.0;
1995	at = 0.0;
1996	}
1997	} // PK
1998
1999	// special treatment for unbound nuclei
2000	if (sn < 0.0) {
2001	probn = 1.0;
2002	probp = 0.0;
2003	proba = 0.0;
2004	probf = 0.0;
2005	goto direct70;
2006	}
2007	if (sbp < 0.0) {
2008	probp = 1.0;
2009	probn = 0.0;
2010	proba = 0.0;
2011	probf = 0.0;
2012	goto direct70;
2013	}
2014
2015	if ((a < 50.e0) \|\| (ee > edyn)) { // no fission if e*> edyn or mass < 50
2016	// G4cout <<"densf = 0.0" << G4endl;
2017	densf = 0.e0;
2018	}
2019
2020	bshell = ecld->ecgnz[in][iz] - ecld->vgsld[in][iz];
2021	defbet = 1.5e0 * (ecld->alpha[in][iz]);
2022
2023	// compound nucleus level density
2024	densniv(a,zprf,ee,0.0e0,&densg,bshell,1.e0,1.e0,&temp,int(fiss->optshp),int(fiss->optcol),defbet);
2025	// assert(isnan(densg) == false);
2026	// assert(isnan(temp) == false);
2027
2028	if ( densg > 0.e0) {
2029	// calculation of the partial decay width
2030	// used for both the time scale and the evaporation decay width
2031	gp = (std::pow(a,(2.0/3.0))/fiss->akap)densp/densg/pistd::pow(pt,2);
2032	gn = (std::pow(a,(2.0/3.0))/fiss->akap)densn/densg/pistd::pow(nt,2);
2033	ga = (std::pow(a,(2.0/3.0))/fiss->akap)densa/densg/pi2.0*std::pow(at,2);
2034	gf = densf/densg/pi/2.0*ft;
2035	// assert(isnan(gf) == false);
2036
2037	// assert(isnan(gp) == false);
2038	// assert(isnan(gn) == false);
2039	// assert(isnan(ga) == false);
2040	// assert(isnan(ft) == false);
2041	// assert(ft != 0);
2042	// assert(isnan(gf) == false);
2043
2044	if(itest == 1) {
2045	G4cout <<"gn,gp,ga,gf " << gn <<","<< gp <<","<< ga <<","<< gf << G4endl;
2046	}
2047	}
2048	else {
2049	if(verboseLevel > 2) {
2050	G4cout <<"direct: densg <= 0.e0 " << a <<","<< zprf <<","<< ee << G4endl;
2051	}
2052	}
2053
2054	gsum = ga + gp + gn;
2055	// assert(isinf(gsum) == false);
2056	// assert(isnan(gsum) == false);
2057	if (gsum > 0.0) {
2058	ts1 = hbar / gsum;
2059	}
2060	else {
2061	ts1 = 1.0e99;
2062	}
2063
2064	if (inum > ilast) { // new event means reset the time scale
2065	tsum = 0;
2066	}
2067
2068	// calculate the relative probabilities for all decay channels
2069	if (densf == 0.0) {
2070	if (densp == 0.0) {
2071	if (densn == 0.0) {
2072	if (densa == 0.0) {
2073	// no reaction is possible
2074	probf = 0.0;
2075	probp = 0.0;
2076	probn = 0.0;
2077	proba = 0.0;
2078	goto direct70;
2079	}
2080
2081	// alpha evaporation is the only open channel
2082	rf = 0.0;
2083	rp = 0.0;
2084	rn = 0.0;
2085	ra = 1.0;
2086	goto direct50;
2087	}
2088
2089	// alpha emission and neutron emission
2090	rf = 0.0;
2091	rp = 0.0;
2092	rn = 1.0;
2093	ra = densa2.0/densnstd::pow((at/nt),2);
2094	goto direct50;
2095	}
2096	// alpha, proton and neutron emission
2097	rf = 0.0;
2098	rp = 1.0;
2099	rn = densn/densp*std::pow((nt/pt),2);
2100	// assert(isnan(rn) == false);
2101	ra = densa2.0/denspstd::pow((at/pt),2);
2102	// assert(isnan(ra) == false);
2103	goto direct50;
2104	}
2105
2106	// here fission has taken place
2107	rf = 1.0;
2108
2109	// cramers and weidenmueller factors for the dynamical hindrances of
2110	// fission
2111	if (fiss->bet <= 1.0e-16) {
2112	cf = 1.0;
2113	wf = 1.0;
2114	}
2115	else if (sbf > 0.0e0) {
2116	cf = cram(fiss->bet,fiss->homega);
2117	// if fission barrier ef=0.d0 then fission is the only possible
2118	// channel. to avoid std::log(0) in function tau
2119	// a.j. 7/28/93
2120	if (ef <= 0.0) {
2121	rp = 0.0;
2122	rn = 0.0;
2123	ra = 0.0;
2124	goto direct50;
2125	}
2126	else {
2127	// transient time tau()
2128	tauc = tau(fiss->bet,fiss->homega,ef,ft);
2129	// assert(isnan(tauc) == false);
2130	}
2131	wfex = (tauc - tsum)/ts1;
2132
2133	if (wfex < 0.0) {
2134	wf = 1.0;
2135	}
2136	else {
2137	wf = std::exp( -wfex);
2138	}
2139	}
2140	else {
2141	cf=1.0;
2142	wf=1.0;
2143	}
2144
2145	if(verboseLevel > 2) {
2146	G4cout <<"tsum,wf,cf " << tsum <<","<< wf <<","<< cf << G4endl;
2147	}
2148
2149	tsum = tsum + ts1;
2150
2151	// change by g.k. and a.h. 5.9.95
2152	tconst = 0.7;
2153	dconst = 12.0/std::sqrt(a);
2154	// assert(isnan(dconst) == false);
2155	nprf = a - zprf;
2156
2157	if (fiss->optshp >= 2) { //then
2158	parite(nprf,&parc);
2159	// assert(isnan(parc) == false);
2160	dconst = dconst*parc;
2161	}
2162	else {
2163	dconst= 0.0;
2164	}
2165	if ((ee <= 17.e0) && (fiss->optles == 1) && (iz >= 90) && (in >= 134)) { //then
2166	// constant changed to 5.0 accord to moretto & vandenbosch a.j. 19.3.96
2167	gngf = std::pow(a,(2.0/3.0))tconst/10.0std::exp((ef-sn+dconst)/tconst);
2168
2169	// if the excitation energy is so low that densn=0 ==> gn = 0
2170	// fission remains the only channel.
2171	// a. j. 10.1.94
2172	if (gn == 0.0) {
2173	rn = 0.0;
2174	rp = 0.0;
2175	ra = 0.0;
2176	}
2177	else {
2178	rn=gngf;
2179	// assert(isnan(rn) == false);
2180	rp=gngf*gp/gn;
2181	// assert(isnan(rp) == false);
2182	ra=gngf*ga/gn;
2183	// assert(isnan(ra) == false);
2184	}
2185	} else {
2186	// assert(isnan(cf) == false);
2187	// assert(isinf(gn) == false);
2188	// assert(isinf(gf) == false);
2189	// assert(isinf(cf) == false);
2190	assert(gn > 0 \|\| (gf != 0 && cf != 0));
2191	rn = gn/(gf*cf);
2192	// G4cout <<"rn = " << G4endl;
2193	// G4cout <<"gn = " << gn << " gf = " << gf << " cf = " << cf << G4endl;
2194	// assert(isnan(rn) == false);
2195	rp = gp/(gf*cf);
2196	// assert(isnan(rp) == false);
2197	ra = ga/(gf*cf);
2198	// assert(isnan(ra) == false);
2199	}
2200	direct50:
2201	// relative decay probabilities
2202	// assert(isnan(ra) == false);
2203	// assert(isnan(rp) == false);
2204	// assert(isnan(rn) == false);
2205	// assert(isnan(rf) == false);
2206
2207	denomi = rp+rn+ra+rf;
2208	// assert(isnan(denomi) == false);
2209	assert(denomi > 0);
2210	// decay probabilities after transient time
2211	probf = rf/denomi;
2212	// assert(isnan(probf) == false);
2213	probp = rp/denomi;
2214	// assert(isnan(probp) == false);
2215	probn = rn/denomi;
2216	// assert(isnan(probn) == false);
2217	proba = ra/denomi;
2218	// assert(isnan(proba) == false);
2219	// assert(isinf(proba) == false);
2220
2221	// new treatment of grange-weidenmueller factor, 5.1.2000, khs !!!
2222
2223	// decay probabilites with transient time included
2224	// assert(isnan(wf) == false);
2225	assert(std::fabs(probf) <= 1.0);
2226	probf = probf * wf;
2227	if(probf == 1.0) {
2228	probp = 0.0;
2229	probn = 0.0;
2230	proba = 0.0;
2231	}
2232	else {
2233	probp = probp * (wf + (1.e0-wf)/(1.e0-probf));
2234	probn = probn * (wf + (1.e0-wf)/(1.e0-probf));
2235	proba = proba * (wf + (1.e0-wf)/(1.e0-probf));
2236	}
2237	direct70:
2238	// assert(isnan(probp) == false);
2239	// assert(isnan(probn) == false);
2240	// assert(isnan(probf) == false);
2241	// assert(isnan(proba) == false);
2242	ptotl = probp+probn+proba+probf;
2243	// assert(isnan(ptotl) == false);
2244
2245	ee = eer;
2246	ilast = inum;
2247
2248	// Return values:
2249	// assert(isnan(proba) == false);
2250	(*probp_par) = probp;
2251	(*probn_par) = probn;
2252	(*proba_par) = proba;
2253	(*probf_par) = probf;
2254	(*ptotl_par) = ptotl;
2255	(*sn_par) = sn;
2256	(*sbp_par) = sbp;
2257	(*sba_par) = sba;
2258	(*ecn_par) = ecn;
2259	(*ecp_par) = ecp;
2260	(*eca_par) = eca;
2261	(*bp_par) = bp;
2262	(*ba_par) = ba;
2263	}
2264
2265	void G4Abla::densniv(G4double a, G4double z, G4double ee, G4double esous, G4double *dens, G4double bshell, G4double bs, G4double bk,
2266	G4double *temp, G4int optshp, G4int optcol, G4double defbet)
2267	{
2268	// 1498 C
2269	// 1499 C INPUT:
2270	// 1500 C A,EE,ESOUS,OPTSHP,BS,BK,BSHELL,DEFBET
2271	// 1501 C
2272	// 1502 C LEVEL DENSITY PARAMETERS
2273	// 1503 C COMMON /ALD/ AV,AS,AK,OPTAFAN
2274	// 1504 C AV,AS,AK - VOLUME,SURFACE,CURVATURE DEPENDENCE OF THE
2275	// 1505 C LEVEL DENSITY PARAMETER
2276	// 1506 C OPTAFAN - 0/1 AF/AN >=1 OR AF/AN ==1
2277	// 1507 C RECOMMENDED IS OPTAFAN = 0
2278	// 1508 C---------------------------------------------------------------------
2279	// 1509 C OUTPUT: DENS,TEMP
2280	// 1510 C
2281	// 1511 C ____________________________________________________________________
2282	// 1512 C /
2283	// 1513 C / PROCEDURE FOR CALCULATING THE STATE DENSITY OF A COMPOUND NUCLEUS
2284	// 1514 C /____________________________________________________________________
2285	// 1515 C
2286	// 1516 INTEGER AFP,IZ,OPTSHP,OPTCOL,J,OPTAFAN
2287	// 1517 REAL*8 A,EE,ESOUS,DENS,E,Y0,Y1,Y2,Y01,Y11,Y21,PA,BS,BK,TEMP
2288	// 1518 C=====INSERTED BY KUDYAEV===============================================
2289	// 1519 COMMON /ALD/ AV,AS,AK,OPTAFAN
2290	// 1520 REAL*8 ECR,ER,DELTAU,Z,DELTPP,PARA,PARZ,FE,HE,ECOR,ECOR1,Pi6
2291	// 1521 REAL*8 BSHELL,DELTA0,AV,AK,AS,PONNIV,PONFE,DEFBET,QR,SIG,FP
2292	// 1522 C=======================================================================
2293	// 1523 C
2294	// 1524 C
2295	// 1525 C-----------------------------------------------------------------------
2296	// 1526 C A MASS NUMBER OF THE DAUGHTER NUCLEUS
2297	// 1527 C EE EXCITATION ENERGY OF THE MOTHER NUCLEUS
2298	// 1528 C ESOUS SEPARATION ENERGY PLUS EFFECTIVE COULOMB BARRIER
2299	// 1529 C DENS STATE DENSITY OF DAUGHTER NUCLEUS AT EE-ESOUS-EC
2300	// 1530 C BSHELL SHELL CORRECTION
2301	// 1531 C TEMP NUCLEAR TEMPERATURE
2302	// 1532 C E LOCAL EXCITATION ENERGY OF THE DAUGHTER NUCLEUS
2303	// 1533 C E1 LOCAL HELP VARIABLE
2304	// 1534 C Y0,Y1,Y2,Y01,Y11,Y21
2305	// 1535 C LOCAL HELP VARIABLES
2306	// 1536 C PA LOCAL STATE-DENSITY PARAMETER
2307	// 1537 C EC KINETIC ENERGY OF EMITTED PARTICLE WITHOUT
2308	// 1538 C COULOMB REPULSION
2309	// 1539 C IDEN FAKTOR FOR SUBSTRACTING KINETIC ENERGY IDEN*TEMP
2310	// 1540 C DELTA0 PAIRING GAP 12 FOR GROUND STATE
2311	// 1541 C 14 FOR SADDLE POINT
2312	// 1542 C EITERA HELP VARIABLE FOR TEMPERATURE ITERATION
2313	// 1543 C-----------------------------------------------------------------------
2314	// 1544 C
2315	// 1545 C
2316	G4double afp = 0.0;
2317	G4double delta0 = 0.0;
2318	G4double deltau = 0.0;
2319	G4double deltpp = 0.0;
2320	G4double e = 0.0;
2321	G4double ecor = 0.0;
2322	G4double ecor1 = 0.0;
2323	G4double ecr = 0.0;
2324	G4double er = 0.0;
2325	G4double fe = 0.0;
2326	G4double fp = 0.0;
2327	G4double he = 0.0;
2328	G4double iz = 0.0;
2329	G4double pa = 0.0;
2330	G4double para = 0.0;
2331	G4double parz = 0.0;
2332	G4double ponfe = 0.0;
2333	G4double ponniv = 0.0;
2334	G4double qr = 0.0;
2335	G4double sig = 0.0;
2336	G4double y01 = 0.0;
2337	G4double y11 = 0.0;
2338	G4double y2 = 0.0;
2339	G4double y21 = 0.0;
2340	G4double y1 = 0.0;
2341	G4double y0 = 0.0;
2342
2343	G4double pi6 = std::pow(3.1415926535,2) / 6.0;
2344	ecr=10.0;
2345	er=28.0;
2346	afp=idnint(a);
2347	iz=idnint(z);
2348
2349	// level density parameter
2350	if((ald->optafan == 1)) {
2351	pa = (ald->av)a + (ald->as)std::pow(a,(2.e0/3.e0)) + (ald->ak)*std::pow(a,(1.e0/3.e0));
2352	}
2353	else {
2354	pa = (ald->av)a + (ald->as)bsstd::pow(a,(2.e0/3.e0)) + (ald->ak)bk*std::pow(a,(1.e0/3.e0));
2355	}
2356
2357	fp = 0.01377937231e0 * std::pow(a,(5.e0/3.e0)) * (1.e0 + defbet/3.e0);
2358
2359	// pairing corrections
2360	if (bs > 1.0) {
2361	delta0 = 14;
2362	}
2363	else {
2364	delta0 = 12;
2365	}
2366
2367	if (esous > 1.0e30) {
2368	(*dens) = 0.0;
2369	(*temp) = 0.0;
2370	goto densniv100;
2371	}
2372
2373	e = ee - esous;
2374
2375	if (e < 0.0) {
2376	(*dens) = 0.0;
2377	(*temp) = 0.0;
2378	goto densniv100;
2379	}
2380
2381	// shell corrections
2382	if (optshp > 0) {
2383	deltau = bshell;
2384	if (optshp == 2) {
2385	deltau = 0.0;
2386	}
2387	if (optshp >= 2) {
2388	// pairing energy shift with condensation energy a.r.j. 10.03.97
2389	// deltpp = -0.25e0* (delta0/std::pow(std::sqrt(a),2)) * pa /pi6 + 2.e0*delta0/std::sqrt(a);
2390	deltpp = -0.25e0* std::pow((delta0/std::sqrt(a)),2) * pa /pi6 + 2.e0*delta0/std::sqrt(a);
2391	// assert(isnan(deltpp) == false);
2392
2393	parite(a,&para);
2394	if (para < 0.0) {
2395	e = e - delta0/std::sqrt(a);
2396	// assert(isnan(e) == false);
2397	} else {
2398	parite(z,&parz);
2399	if (parz > 0.e0) {
2400	e = e - 2.0*delta0/std::sqrt(a);
2401	// assert(isnan(e) == false);
2402	} else {
2403	e = e;
2404	// assert(isnan(e) == false);
2405	}
2406	}
2407	} else {
2408	deltpp = 0.0;
2409	}
2410	} else {
2411	deltau = 0.0;
2412	deltpp = 0.0;
2413	}
2414	if (e < 0.0) {
2415	e = 0.0;
2416	(*temp) = 0.0;
2417	}
2418
2419	// washing out is made stronger ! g.k. 3.7.96
2420	ponfe = -2.5pae*std::pow(a,(-4.0/3.0));
2421
2422	if (ponfe < -700.0) {
2423	ponfe = -700.0;
2424	}
2425	fe = 1.0 - std::exp(ponfe);
2426	if (e < ecr) {
2427	// priv. comm. k.-h. schmidt
2428	he = 1.0 - std::pow((1.0 - e/ecr),2);
2429	}
2430	else {
2431	he = 1.0;
2432	}
2433
2434	// Excitation energy corrected for pairing and shell effects
2435	// washing out with excitation energy is included.
2436	ecor = e + deltaufe + deltpphe;
2437
2438	if (ecor <= 0.1) {
2439	ecor = 0.1;
2440	}
2441
2442	// statt 170.d0 a.r.j. 8.11.97
2443
2444	// iterative procedure according to grossjean and feldmeier
2445	// to avoid the singularity e = 0
2446	if (ee < 5.0) {
2447	y1 = std::sqrt(pa*ecor);
2448	// assert(isnan(y1) == false);
2449	for(int j = 0; j < 5; j++) {
2450	y2 = paecor(1.e0-std::exp(-y1));
2451	// assert(isnan(y2) == false);
2452	y1 = std::sqrt(y2);
2453	// assert(isnan(y1) == false);
2454	}
2455
2456	y0 = pa/y1;
2457	// assert(isnan(y0) == false);
2458	assert(y0 != 0.0);
2459	(*temp)=1.0/y0;
2460	(dens) = std::exp(y0ecor)/ (std::pow((std::pow(ecor,3)y0),0.5)std::pow((1.0-0.5y0ecorstd::exp(-y1)),0.5)) std::exp(y1)(1.0-std::exp(-y1))0.1477045;
2461	if (ecor < 1.0) {
2462	ecor1=1.0;
2463	y11 = std::sqrt(pa*ecor1);
2464	// assert(isnan(y11) == false);
2465	for(int j = 0; j < 7; j++) {
2466	y21 = paecor1(1.0-std::exp(-y11));
2467	// assert(isnan(21) == false);
2468	y11 = std::sqrt(y21);
2469	// assert(isnan(y11) == false);
2470	}
2471
2472	y01 = pa/y11;
2473	// assert(isnan(y01) == false);
2474	(dens) = (dens)*std::pow((y01/y0),1.5);
2475	(temp) = (temp)*std::pow((y01/y0),1.5);
2476	}
2477	}
2478	else {
2479	ponniv = 2.0std::sqrt(paecor);
2480	// assert(isnan(ponniv) == false);
2481	if (ponniv > 700.0) {
2482	ponniv = 700.0;
2483	}
2484
2485	// fermi gas state density
2486	(dens) = std::pow(pa,(-0.25e0))std::pow(ecor,(-1.25e0))std::exp(ponniv) 0.1477045e0;
2487	// assert(isnan(std::sqrt(ecor/pa)) == false);
2488	(*temp) = std::sqrt(ecor/pa);
2489	}
2490	densniv100:
2491
2492	// spin cutoff parameter
2493	sig = fp * (*temp);
2494
2495	// collective enhancement
2496	if (optcol == 1) {
2497	qrot(z,a,defbet,sig,ecor,&qr);
2498	}
2499	else {
2500	qr = 1.0;
2501	}
2502
2503	(dens) = (dens) * qr;
2504	}
2505
2506
2507	G4double G4Abla::bfms67(G4double zms, G4double ams)
2508	{
2509	// This subroutine calculates the fission barriers
2510	// of the liquid-drop model of Myers and Swiatecki (1967).
2511	// Analytic parameterization of Dahlinger 1982
2512	// replaces tables. Barrier heights from Myers and Swiatecki !!!
2513
2514	G4double nms = 0.0, ims = 0.0, ksims = 0.0, xms = 0.0, ums = 0.0;
2515
2516	nms = ams - zms;
2517	ims = (nms-zms)/ams;
2518	ksims= 50.15e0 * (1.- 1.78 * std::pow(ims,2));
2519	xms = std::pow(zms,2) / (ams * ksims);
2520	ums = 0.368e0-5.057e0xms+8.93e0std::pow(xms,2)-8.71*std::pow(xms,3);
2521	return(0.7322e0std::pow(zms,2)/std::pow(ams,(0.333333e0))std::pow(10.e0,ums));
2522	}
2523
2524	void G4Abla::lpoly(G4double x, G4int n, G4double pl[])
2525	{
2526	// THIS SUBROUTINE CALCULATES THE ORDINARY LEGENDRE POLYNOMIALS OF
2527	// ORDER 0 TO N-1 OF ARGUMENT X AND STORES THEM IN THE VECTOR PL.
2528	// THEY ARE CALCULATED BY RECURSION RELATION FROM THE FIRST TWO
2529	// POLYNOMIALS.
2530	// WRITTEN BY A.J.SIERK LANL T-9 FEBRUARY, 1984
2531
2532	// NOTE: PL AND X MUST BE DOUBLE PRECISION ON 32-BIT COMPUTERS!
2533
2534	pl[0] = 1.0;
2535	pl[1] = x;
2536
2537	for(int i = 2; i < n; i++) {
2538	pl[i] = ((2double(i+1) - 3.0)xpl[i-1] - (double(i+1) - 2.0)pl[i-2])/(double(i+1)-1.0);
2539	}
2540	}
2541
2542	G4double G4Abla::eflmac(G4int ia, G4int iz, G4int flag, G4int optshp)
2543	{
2544	// CHANGED TO CALCULATE TOTAL BINDING ENERGY INSTEAD OF MASS EXCESS.
2545	// SWITCH FOR PAIRING INCLUDED AS WELL.
2546	// BINDING = EFLMAC(IA,IZ,0,OPTSHP)
2547	// FORTRAN TRANSCRIPT OF /U/GREWE/LANG/EEX/FRLDM.C
2548	// A.J. 15.07.96
2549
2550	// this function will calculate the liquid-drop nuclear mass for spheri
2551	// configuration according to the preprint NUCLEAR GROUND-STATE
2552	// MASSES and DEFORMATIONS by P. M"oller et al. from August 16, 1993 p.
2553	// All constants are taken from this publication for consistency.
2554
2555	// Parameters:
2556	// a: nuclear mass number
2557	// z: nuclear charge
2558	// flag: 0 - return mass excess
2559	// otherwise - return pairing (= -1/2 dpn + 1/2 (Dp + Dn))
2560
2561	G4double eflmacResult = 0.0;
2562
2563	G4int in = 0;
2564	G4double z = 0.0, n = 0.0, a = 0.0, av = 0.0, as = 0.0;
2565	G4double a0 = 0.0, c1 = 0.0, c4 = 0.0, b1 = 0.0, b3 = 0.0;
2566	G4double f = 0.0, ca = 0.0, w = 0.0, dp = 0.0, dn = 0.0, dpn = 0.0, efl = 0.0;
2567	G4double rmac = 0.0, bs = 0.0, h = 0.0, r0 = 0.0, kf = 0.0, ks = 0.0;
2568	G4double kv = 0.0, rp = 0.0, ay = 0.0, aden = 0.0, x0 = 0.0, y0 = 0.0;
2569	G4double mh = 0.0, mn = 0.0, esq = 0.0, ael = 0.0, i = 0.0;
2570	G4double pi = 3.141592653589793238e0;
2571
2572	// fundamental constants
2573	// hydrogen-atom mass excess
2574	mh = 7.289034;
2575
2576	// neutron mass excess
2577	mn = 8.071431;
2578
2579	// electronic charge squared
2580	esq = 1.4399764;
2581
2582	// constants from considerations other than nucl. masses
2583	// electronic binding
2584	ael = 1.433e-5;
2585
2586	// proton rms radius
2587	rp = 0.8;
2588
2589	// nuclear radius constant
2590	r0 = 1.16;
2591
2592	// range of yukawa-plus-expon. potential
2593	ay = 0.68;
2594
2595	// range of yukawa function used to generate
2596	// nuclear charge distribution
2597	aden= 0.70;
2598
2599	// constants from considering odd-even mass differences
2600	// average pairing gap
2601	rmac= 4.80;
2602
2603	// neutron-proton interaction
2604	h = 6.6;
2605
2606	// wigner constant
2607	w = 30.0;
2608
2609	// adjusted parameters
2610	// volume energy
2611	av = 16.00126;
2612
2613	// volume asymmetry
2614	kv = 1.92240;
2615
2616	// surface energy
2617	as = 21.18466;
2618
2619	// surface asymmetry
2620	ks = 2.345;
2621	// a^0 constant
2622	a0 = 2.615;
2623
2624	// charge asymmetry
2625	ca = 0.10289;
2626
2627	// we will account for deformation by using the microscopic
2628	// corrections tabulated from p. 68ff */
2629	bs = 1.0;
2630
2631	z = double(iz);
2632	a = double(ia);
2633	in = ia - iz;
2634	n = double(in);
2635	dn = rmac*bs/std::pow(n,(1.0/3.0));
2636	dp = rmac*bs/std::pow(z,(1.0/3.0));
2637	dpn = h/bs/std::pow(a,(2.0/3.0));
2638	// assert(isnan(dpn) == false);
2639
2640	c1 = 3.0/5.0*esq/r0;
2641	// assert(isnan(c1) == false);
2642	// assert(isinf(c1) == false);
2643
2644	c4 = 5.0/4.0std::pow((3.0/(2.0pi)),(2.0/3.0)) * c1;
2645	// assert(isnan(c4) == false);
2646	// assert(isinf(c4) == false);
2647
2648	// assert(isnan(pi) == false);
2649	// assert(isnan(z) == false);
2650	// assert(isnan(a) == false);
2651	// assert(isnan(r0) == false);
2652	kf = std::pow((9.0piz/(4.0*a)),(1.0/3.0))/r0;
2653	// assert(isnan(kf) == false);
2654	// assert(isinf(kf) == false);
2655
2656	f = -1.0/8.0rprpesq/std::pow(r0,3) (145.0/48.0 - 327.0/2880.0std::pow(kf,2) std::pow(rp,2) + 1527.0/1209600.0std::pow(kf,4) std::pow(rp,4));
2657	i = (n-z)/a;
2658
2659	x0 = r0 * std::pow(a,(1.0/3.0)) / ay;
2660	y0 = r0 * std::pow(a,(1.0/3.0)) / aden;
2661
2662	b1 = 1.0 - 3.0/(std::pow(x0,2)) + (1.0 + x0) * (2.0 + 3.0/x0 + 3.0/std::pow(x0,2)) * std::exp(-2.0*x0);
2663
2664	b3 = 1.0 - 5.0/std::pow(y0,2) * (1.0 - 15.0/(8.0y0) + 21.0/(8.0 std::pow(y0,3))
2665	- 3.0/4.0 * (1.0 + 9.0/(2.0*y0) + 7.0/std::pow(y0,2)
2666	+ 7.0/(2.0 * std::pow(y0,3))) * std::exp(-2.0*y0));
2667
2668	// now calulation of total binding energy a.j. 16.7.96
2669
2670	efl = -1.0 * av(1.0 - kvii)a + as(1.0 - ksii)b1 * std::pow(a,(2.0/3.0)) + a0
2671	+ c1zzb3/std::pow(a,(1.0/3.0)) - c4std::pow(z,(4.0/3.0))/std::pow(a,(1.e0/3.e0))
2672	+ fstd::pow(z,2)/a -ca(n-z) - ael * std::pow(z,(2.39e0));
2673
2674	if ((in == iz) && (mod(in,2) == 1) && (mod(iz,2) == 1)) {
2675	// n and z odd and equal
2676	efl = efl + w*(utilabs(i)+1.e0/a);
2677	}
2678	else {
2679	efl= efl + w* utilabs(i);
2680	}
2681
2682	// pairing is made optional
2683	if (optshp >= 2) {
2684	// average pairing
2685	if ((mod(in,2) == 1) && (mod(iz,2) == 1)) {
2686	efl = efl - dpn;
2687	}
2688	if (mod(in,2) == 1) {
2689	efl = efl + dn;
2690	}
2691	if (mod(iz,2) == 1) {
2692	efl = efl + dp;
2693	}
2694	// end if for pairing term
2695	}
2696
2697	if (flag != 0) {
2698	eflmacResult = (0.5(dn + dp) - 0.5dpn);
2699	}
2700	else {
2701	eflmacResult = efl;
2702	}
2703
2704	return eflmacResult;
2705	}
2706
2707	void G4Abla::appariem(G4double a, G4double z, G4double *del)
2708	{
2709	// CALCUL DE LA CORRECTION, DUE A L'APPARIEMENT, DE L'ENERGIE DE
2710	// LIAISON D'UN NOYAU
2711	// PROCEDURE FOR CALCULATING THE PAIRING CORRECTION TO THE BINDING
2712	// ENERGY OF A SPECIFIC NUCLEUS
2713
2714	double para = 0.0, parz = 0.0;
2715	// A MASS NUMBER
2716	// Z NUCLEAR CHARGE
2717	// PARA HELP VARIABLE FOR PARITY OF A
2718	// PARZ HELP VARIABLE FOR PARITY OF Z
2719	// DEL PAIRING CORRECTION
2720
2721	parite(a, &para);
2722
2723	if (para < 0.0) {
2724	(*del) = 0.0;
2725	}
2726	else {
2727	parite(z, &parz);
2728	if (parz > 0.0) {
2729	// assert(isnan(std::sqrt(a)) == false);
2730	(*del) = -12.0/std::sqrt(a);
2731	}
2732	else {
2733	// assert(isnan(std::sqrt(a)) == false);
2734	(*del) = 12.0/std::sqrt(a);
2735	}
2736	}
2737	}
2738
2739	void G4Abla::parite(G4double n, G4double *par)
2740	{
2741	// CALCUL DE LA PARITE DU NOMBRE N
2742	//
2743	// PROCEDURE FOR CALCULATING THE PARITY OF THE NUMBER N.
2744	// RETURNS -1 IF N IS ODD AND +1 IF N IS EVEN
2745
2746	G4double n1 = 0.0, n2 = 0.0, n3 = 0.0;
2747
2748	// N NUMBER TO BE TESTED
2749	// N1,N2 HELP VARIABLES
2750	// PAR HELP VARIABLE FOR PARITY OF N
2751
2752	n3 = double(idnint(n));
2753	n1 = n3/2.0;
2754	n2 = n1 - dint(n1);
2755
2756	if (n2 > 0.0) {
2757	(*par) = -1.0;
2758	}
2759	else {
2760	(*par) = 1.0;
2761	}
2762	}
2763
2764	G4double G4Abla::tau(G4double bet, G4double homega, G4double ef, G4double t)
2765	{
2766	// INPUT : BET, HOMEGA, EF, T
2767	// OUTPUT: TAU - RISE TIME IN WHICH THE FISSION WIDTH HAS REACHED
2768	// 90 PERCENT OF ITS FINAL VALUE
2769	//
2770	// BETA - NUCLEAR VISCOSITY
2771	// HOMEGA - CURVATURE OF POTENTIAL
2772	// EF - FISSION BARRIER
2773	// T - NUCLEAR TEMPERATURE
2774
2775	G4double tauResult = 0.0;
2776
2777	G4double tlim = 8.e0 * ef;
2778	if (t > tlim) {
2779	t = tlim;
2780	}
2781
2782	// modified bj and khs 6.1.2000:
2783	if (bet/(std::sqrt(2.0)10.0(homega/6.582122)) <= 1.0) {
2784	tauResult = std::log(10.0ef/t)/(bet1.0e21);
2785	// assert(isnan(tauResult) == false);
2786	}
2787	else {
2788	tauResult = std::log(10.0ef/t)/ (2.0std::pow((10.0homega/6.582122),2))(bet*1.0e-21);
2789	// assert(isnan(tauResult) == false);
2790	} //end if
2791
2792	return tauResult;
2793	}
2794
2795	G4double G4Abla::cram(G4double bet, G4double homega)
2796	{
2797	// INPUT : BET, HOMEGA NUCLEAR VISCOSITY + CURVATURE OF POTENTIAL
2798	// OUTPUT: KRAMERS FAKTOR - REDUCTION OF THE FISSION PROBABILITY
2799	// INDEPENDENT OF EXCITATION ENERGY
2800
2801	G4double rel = bet/(20.0*homega/6.582122);
2802	G4double cramResult = std::sqrt(1.0 + std::pow(rel,2)) - rel;
2803	// limitation introduced 6.1.2000 by khs
2804
2805	if (cramResult > 1.0) {
2806	cramResult = 1.0;
2807	}
2808
2809	// assert(isnan(cramResult) == false);
2810	return cramResult;
2811	}
2812
2813	G4double G4Abla::bipol(int iflag, G4double y)
2814	{
2815	// CALCULATION OF THE SURFACE BS OR CURVATURE BK OF A NUCLEUS
2816	// RELATIVE TO THE SPHERICAL CONFIGURATION
2817	// BASED ON MYERS, DROPLET MODEL FOR ARBITRARY SHAPES
2818
2819	// INPUT: IFLAG - 0/1 BK/BS CALCULATION
2820	// Y - (1 - X) COMPLEMENT OF THE FISSILITY
2821
2822	// LINEAR INTERPOLATION OF BS BK TABLE
2823
2824	int i = 0;
2825
2826	G4double bipolResult = 0.0;
2827
2828	const int bsbkSize = 54;
2829
2830	G4double bk[bsbkSize] = {0.0, 1.00000,1.00087,1.00352,1.00799,1.01433,1.02265,1.03306,
2831	1.04576,1.06099,1.07910,1.10056,1.12603,1.15651,1.19348,
2832	1.23915,1.29590,1.35951,1.41013,1.44103,1.46026,1.47339,
2833	1.48308,1.49068,1.49692,1.50226,1.50694,1.51114,1.51502,
2834	1.51864,1.52208,1.52539,1.52861,1.53177,1.53490,1.53803,
2835	1.54117,1.54473,1.54762,1.55096,1.55440,1.55798,1.56173,
2836	1.56567,1.56980,1.57413,1.57860,1.58301,1.58688,1.58688,
2837	1.58688,1.58740,1.58740, 0.0}; //Zeroes at bk[0], and at the end added by PK
2838
2839	G4double bs[bsbkSize] = {0.0, 1.00000,1.00086,1.00338,1.00750,1.01319,
2840	1.02044,1.02927,1.03974,
2841	1.05195,1.06604,1.08224,1.10085,1.12229,1.14717,1.17623,1.20963,
2842	1.24296,1.26532,1.27619,1.28126,1.28362,1.28458,1.28477,1.28450,
2843	1.28394,1.28320,1.28235,1.28141,1.28042,1.27941,1.27837,1.27732,
2844	1.27627,1.27522,1.27418,1.27314,1.27210,1.27108,1.27006,1.26906,
2845	1.26806,1.26707,1.26610,1.26514,1.26418,1.26325,1.26233,1.26147,
2846	1.26147,1.26147,1.25992,1.25992, 0.0};
2847
2848	i = idint(y/(2.0e-02)) + 1;
2849	assert(i >= 1);
2850
2851	if(i >= bsbkSize) {
2852	if(verboseLevel > 2) {
2853	G4cout <<"G4Abla error: index i = " << i << " is greater than array size permits." << G4endl;
2854	}
2855	bipolResult = 0.0;
2856	}
2857	else {
2858	if (iflag == 1) {
2859	bipolResult = bs[i] + (bs[i+1] - bs[i])/2.0e-02 * (y - 2.0e-02*(i - 1));
2860	}
2861	else {
2862	bipolResult = bk[i] + (bk[i+1] - bk[i])/2.0e-02 * (y - 2.0e-02*(i - 1));
2863	}
2864	}
2865
2866	// assert(isnan(bipolResult) == false);
2867	return bipolResult;
2868	}
2869
2870	void G4Abla::barfit(G4int iz, G4int ia, G4int il, G4double sbfis, G4double segs, G4double *selmax)
2871	{
2872	// 2223 C VERSION FOR 32BIT COMPUTER
2873	// 2224 C THIS SUBROUTINE RETURNS THE BARRIER HEIGHT BFIS, THE
2874	// 2225 C GROUND-STATE ENERGY SEGS, IN MEV, AND THE ANGULAR MOMENTUM
2875	// 2226 C AT WHICH THE FISSION BARRIER DISAPPEARS, LMAX, IN UNITS OF
2876	// 2227 C H-BAR, WHEN CALLED WITH INTEGER AGUMENTS IZ, THE ATOMIC
2877	// 2228 C NUMBER, IA, THE ATOMIC MASS NUMBER, AND IL, THE ANGULAR
2878	// 2229 C MOMENTUM IN UNITS OF H-BAR. (PLANCK'S CONSTANT DIVIDED BY
2879	// 2230 C 2*PI).
2880	// 2231 C
2881	// 2232 C THE FISSION BARRIER FO IL = 0 IS CALCULATED FROM A 7TH
2882	// 2233 C ORDER FIT IN TWO VARIABLES TO 638 CALCULATED FISSION
2883	// 2234 C BARRIERS FOR Z VALUES FROM 20 TO 110. THESE 638 BARRIERS ARE
2884	// 2235 C FIT WITH AN RMS DEVIATION OF 0.10 MEV BY THIS 49-PARAMETER
2885	// 2236 C FUNCTION.
2886	// 2237 C IF BARFIT IS CALLED WITH (IZ,IA) VALUES OUTSIDE THE RANGE OF
2887	// 2238 C THE BARRIER HEIGHT IS SET TO 0.0, AND A MESSAGE IS PRINTED
2888	// 2239 C ON THE DEFAULT OUTPUT FILE.
2889	// 2240 C
2890	// 2241 C FOR IL VALUES NOT EQUAL TO ZERO, THE VALUES OF L AT WHICH
2891	// 2242 C THE BARRIER IS 80% AND 20% OF THE L=0 VALUE ARE RESPECTIVELY
2892	// 2243 C FIT TO 20-PARAMETER FUNCTIONS OF Z AND A, OVER A MORE
2893	// 2244 C RESTRICTED RANGE OF A VALUES, THAN IS THE CASE FOR L = 0.
2894	// 2245 C THE VALUE OF L WHERE THE BARRIER DISAPPEARS, LMAX IS FIT TO
2895	// 2246 C A 24-PARAMETER FUNCTION OF Z AND A, WITH THE SAME RANGE OF
2896	// 2247 C Z AND A VALUES AS L-80 AND L-20.
2897	// 2248 C ONCE AGAIN, IF AN (IZ,IA) PAIR IS OUTSIDE OF THE RANGE OF
2898	// 2249 C VALIDITY OF THE FIT, THE BARRIER VALUE IS SET TO 0.0 AND A
2899	// 2250 C MESSAGE IS PRINTED. THESE THREE VALUES (BFIS(L=0),L-80, AND
2900	// 2251 C L-20) AND THE CONSTRINTS OF BFIS = 0 AND D(BFIS)/DL = 0 AT
2901	// 2252 C L = LMAX AND L=0 LEAD TO A FIFTH-ORDER FIT TO BFIS(L) FOR
2902	// 2253 C L>L-20. THE FIRST THREE CONSTRAINTS LEAD TO A THIRD-ORDER FIT
2903	// 2254 C FOR THE REGION L < L-20.
2904	// 2255 C
2905	// 2256 C THE GROUND STATE ENERGIES ARE CALCULATED FROM A
2906	// 2257 C 120-PARAMETER FIT IN Z, A, AND L TO 214 GROUND-STATE ENERGIES
2907	// 2258 C FOR 36 DIFFERENT Z AND A VALUES.
2908	// 2259 C (THE RANGE OF Z AND A IS THE SAME AS FOR L-80, L-20, AND
2909	// 2260 C L-MAX)
2910	// 2261 C
2911	// 2262 C THE CALCULATED BARRIERS FROM WHICH THE FITS WERE MADE WERE
2912	// 2263 C CALCULATED IN 1983-1984 BY A. J. SIERK OF LOS ALAMOS
2913	// 2264 C NATIONAL LABORATORY GROUP T-9, USING YUKAWA-PLUS-EXPONENTIAL
2914	// 2265 C G4DOUBLE FOLDED NUCLEAR ENERGY, EXACT COULOMB DIFFUSENESS
2915	// 2266 C CORRECTIONS, AND DIFFUSE-MATTER MOMENTS OF INERTIA.
2916	// 2267 C THE PARAMETERS OF THE MODEL R-0 = 1.16 FM, AS 21.13 MEV,
2917	// 2268 C KAPPA-S = 2.3, A = 0.68 FM.
2918	// 2269 C THE DIFFUSENESS OF THE MATTER AND CHARGE DISTRIBUTIONS USED
2919	// 2270 C CORRESPONDS TO A SURFACE DIFFUSENESS PARAMETER (DEFINED BY
2920	// 2271 C MYERS) OF 0.99 FM. THE CALCULATED BARRIERS FOR L = 0 ARE
2921	// 2272 C ACCURATE TO A LITTLE LESS THAN 0.1 MEV; THE OUTPUT FROM
2922	// 2273 C THIS SUBROUTINE IS A LITTLE LESS ACCURATE. WORST ERRORS MAY BE
2923	// 2274 C AS LARGE AS 0.5 MEV; CHARACTERISTIC UNCERTAINY IS IN THE RANGE
2924	// 2275 C OF 0.1-0.2 MEV. THE RMS DEVIATION OF THE GROUND-STATE FIT
2925	// 2276 C FROM THE 214 INPUT VALUES IS 0.20 MEV. THE MAXIMUM ERROR
2926	// 2277 C OCCURS FOR LIGHT NUCLEI IN THE REGION WHERE THE GROUND STATE
2927	// 2278 C IS PROLATE, AND MAY BE GREATER THAN 1.0 MEV FOR VERY NEUTRON
2928	// 2279 C DEFICIENT NUCLEI, WITH L NEAR LMAX. FOR MOST NUCLEI LIKELY TO
2929	// 2280 C BE ENCOUNTERED IN REAL EXPERIMENTS, THE MAXIMUM ERROR IS
2930	// 2281 C CLOSER TO 0.5 MEV, AGAIN FOR LIGHT NUCLEI AND L NEAR LMAX.
2931	// 2282 C
2932	// 2283 C WRITTEN BY A. J. SIERK, LANL T-9
2933	// 2284 C VERSION 1.0 FEBRUARY, 1984
2934	// 2285 C
2935	// 2286 C THE FOLLOWING IS NECESSARY FOR 32-BIT MACHINES LIKE DEC VAX,
2936	// 2287 C IBM, ETC
2937
2938	G4double pa[7],pz[7],pl[10];
2939	for(G4int init_i = 0; init_i < 7; init_i++) {
2940	pa[init_i] = 0.0;
2941	pz[init_i] = 0.0;
2942	}
2943	for(G4int init_i = 0; init_i < 10; init_i++) {
2944	pl[init_i] = 0.0;
2945	}
2946
2947	G4double a = 0.0, z = 0.0, amin = 0.0, amax = 0.0, amin2 = 0.0;
2948	G4double amax2 = 0.0, aa = 0.0, zz = 0.0, bfis = 0.0;
2949	G4double bfis0 = 0.0, ell = 0.0, el = 0.0, egs = 0.0, el80 = 0.0, el20 = 0.0;
2950	G4double elmax = 0.0, sel80 = 0.0, sel20 = 0.0, x = 0.0, y = 0.0, q = 0.0, qa = 0.0, qb = 0.0;
2951	G4double aj = 0.0, ak = 0.0, a1 = 0.0, a2 = 0.0;
2952
2953	G4int i = 0, j = 0, k = 0, m = 0;
2954	G4int l = 0;
2955
2956	G4double emncof[4][5] = {{-9.01100e+2,-1.40818e+3, 2.77000e+3,-7.06695e+2, 8.89867e+2},
2957	{1.35355e+4,-2.03847e+4, 1.09384e+4,-4.86297e+3,-6.18603e+2},
2958	{-3.26367e+3, 1.62447e+3, 1.36856e+3, 1.31731e+3, 1.53372e+2},
2959	{7.48863e+3,-1.21581e+4, 5.50281e+3,-1.33630e+3, 5.05367e-2}};
2960
2961	G4double elmcof[4][5] = {{1.84542e+3,-5.64002e+3, 5.66730e+3,-3.15150e+3, 9.54160e+2},
2962	{-2.24577e+3, 8.56133e+3,-9.67348e+3, 5.81744e+3,-1.86997e+3},
2963	{2.79772e+3,-8.73073e+3, 9.19706e+3,-4.91900e+3, 1.37283e+3},
2964	{-3.01866e+1, 1.41161e+3,-2.85919e+3, 2.13016e+3,-6.49072e+2}};
2965
2966	G4double emxcof[4][6] = {{9.43596e4,-2.241997e5,2.223237e5,-1.324408e5,4.68922e4,-8.83568e3},
2967	{-1.655827e5,4.062365e5,-4.236128e5,2.66837e5,-9.93242e4,1.90644e4},
2968	{1.705447e5,-4.032e5,3.970312e5,-2.313704e5,7.81147e4,-1.322775e4},
2969	{-9.274555e4,2.278093e5,-2.422225e5,1.55431e5,-5.78742e4,9.97505e3}};
2970
2971	G4double elzcof[7][7] = {{5.11819909e+5,-1.30303186e+6, 1.90119870e+6,-1.20628242e+6, 5.68208488e+5, 5.48346483e+4,-2.45883052e+4},
2972	{-1.13269453e+6, 2.97764590e+6,-4.54326326e+6, 3.00464870e+6, -1.44989274e+6,-1.02026610e+5, 6.27959815e+4},
2973	{1.37543304e+6,-3.65808988e+6, 5.47798999e+6,-3.78109283e+6, 1.84131765e+6, 1.53669695e+4,-6.96817834e+4},
2974	{-8.56559835e+5, 2.48872266e+6,-4.07349128e+6, 3.12835899e+6, -1.62394090e+6, 1.19797378e+5, 4.25737058e+4},
2975	{3.28723311e+5,-1.09892175e+6, 2.03997269e+6,-1.77185718e+6, 9.96051545e+5,-1.53305699e+5,-1.12982954e+4},
2976	{4.15850238e+4, 7.29653408e+4,-4.93776346e+5, 6.01254680e+5, -4.01308292e+5, 9.65968391e+4,-3.49596027e+3},
2977	{-1.82751044e+5, 3.91386300e+5,-3.03639248e+5, 1.15782417e+5, -4.24399280e+3,-6.11477247e+3, 3.66982647e+2}};
2978
2979	const G4int sizex = 5;
2980	const G4int sizey = 6;
2981	const G4int sizez = 4;
2982
2983	G4double egscof[sizey][sizey][sizez];
2984
2985	G4double egs1[sizey][sizex] = {{1.927813e5, 7.666859e5, 6.628436e5, 1.586504e5,-7.786476e3},
2986	{-4.499687e5,-1.784644e6,-1.546968e6,-4.020658e5,-3.929522e3},
2987	{4.667741e5, 1.849838e6, 1.641313e6, 5.229787e5, 5.928137e4},
2988	{-3.017927e5,-1.206483e6,-1.124685e6,-4.478641e5,-8.682323e4},
2989	{1.226517e5, 5.015667e5, 5.032605e5, 2.404477e5, 5.603301e4},
2990	{-1.752824e4,-7.411621e4,-7.989019e4,-4.175486e4,-1.024194e4}};
2991
2992	G4double egs2[sizey][sizex] = {{-6.459162e5,-2.903581e6,-3.048551e6,-1.004411e6,-6.558220e4},
2993	{1.469853e6, 6.564615e6, 6.843078e6, 2.280839e6, 1.802023e5},
2994	{-1.435116e6,-6.322470e6,-6.531834e6,-2.298744e6,-2.639612e5},
2995	{8.665296e5, 3.769159e6, 3.899685e6, 1.520520e6, 2.498728e5},
2996	{-3.302885e5,-1.429313e6,-1.512075e6,-6.744828e5,-1.398771e5},
2997	{4.958167e4, 2.178202e5, 2.400617e5, 1.167815e5, 2.663901e4}};
2998
2999	G4double egs3[sizey][sizex] = {{3.117030e5, 1.195474e6, 9.036289e5, 6.876190e4,-6.814556e4},
3000	{-7.394913e5,-2.826468e6,-2.152757e6,-2.459553e5, 1.101414e5},
3001	{7.918994e5, 3.030439e6, 2.412611e6, 5.228065e5, 8.542465e3},
3002	{-5.421004e5,-2.102672e6,-1.813959e6,-6.251700e5,-1.184348e5},
3003	{2.370771e5, 9.459043e5, 9.026235e5, 4.116799e5, 1.001348e5},
3004	{-4.227664e4,-1.738756e5,-1.795906e5,-9.292141e4,-2.397528e4}};
3005
3006	G4double egs4[sizey][sizex] = {{-1.072763e5,-5.973532e5,-6.151814e5, 7.371898e4, 1.255490e5},
3007	{2.298769e5, 1.265001e6, 1.252798e6,-2.306276e5,-2.845824e5},
3008	{-2.093664e5,-1.100874e6,-1.009313e6, 2.705945e5, 2.506562e5},
3009	{1.274613e5, 6.190307e5, 5.262822e5,-1.336039e5,-1.115865e5},
3010	{-5.715764e4,-2.560989e5,-2.228781e5,-3.222789e3, 1.575670e4},
3011	{1.189447e4, 5.161815e4, 4.870290e4, 1.266808e4, 2.069603e3}};
3012
3013	for(i = 0; i < sizey; i++) {
3014	for(j = 0; j < sizex; j++) {
3015	// egscof[i][j][0] = egs1[i][j];
3016	// egscof[i][j][1] = egs2[i][j];
3017	// egscof[i][j][2] = egs3[i][j];
3018	// egscof[i][j][3] = egs4[i][j];
3019	egscof[i][j][0] = egs1[i][j];
3020	egscof[i][j][1] = egs2[i][j];
3021	egscof[i][j][2] = egs3[i][j];
3022	egscof[i][j][3] = egs4[i][j];
3023	}
3024	}
3025
3026	// the program starts here
3027	if (iz < 19 \|\| iz > 111) {
3028	goto barfit900;
3029	}
3030
3031	if(iz > 102 && il > 0) {
3032	goto barfit902;
3033	}
3034
3035	z=double(iz);
3036	a=double(ia);
3037	el=double(il);
3038	amin= 1.2e0z + 0.01e0z*z;
3039	amax= 5.8e0z - 0.024e0z*z;
3040
3041	if(a < amin \|\| a > amax) {
3042	goto barfit910;
3043	}
3044
3045	// angul.mom.zero barrier
3046	aa=2.5e-3*a;
3047	zz=1.0e-2*z;
3048	ell=1.0e-2*el;
3049	bfis0 = 0.0;
3050	lpoly(zz,7,pz);
3051	lpoly(aa,7,pa);
3052
3053	for(i = 0; i < 7; i++) { //do 10 i=1,7
3054	for(j = 0; j < 7; j++) { //do 10 j=1,7
3055	bfis0=bfis0+elzcof[j][i]pz[i]pa[j];
3056	//bfis0=bfis0+elzcof[i][j]pz[j]pa[i];
3057	}
3058	}
3059
3060	bfis=bfis0;
3061	// assert(isnan(bfis) == false);
3062
3063	(*sbfis)=bfis;
3064	egs=0.0;
3065	(*segs)=egs;
3066
3067	// values of l at which the barrier
3068	// is 20%(el20) and 80%(el80) of l=0 value
3069	amin2 = 1.4e0z + 0.009e0z*z;
3070	amax2 = 20.e0 + 3.0e0*z;
3071
3072	if((a < amin2-5.e0 \|\| a > amax2+10.e0) && il > 0) {
3073	goto barfit920;
3074	}
3075
3076	lpoly(zz,5,pz);
3077	lpoly(aa,4,pa);
3078	el80=0.0;
3079	el20=0.0;
3080	elmax=0.0;
3081
3082	for(i = 0; i < 4; i++) {
3083	for(j = 0; j < 5; j++) {
3084	// el80 = el80 + elmcof[j][i]pz[j]pa[i];
3085	// el20 = el20 + emncof[j][i]pz[j]pa[i];
3086	el80 = el80 + elmcof[i][j]pz[j]pa[i];
3087	el20 = el20 + emncof[i][j]pz[j]pa[i];
3088	}
3089	}
3090
3091	sel80 = el80;
3092	sel20 = el20;
3093
3094	// value of l (elmax) where barrier disapp.
3095	lpoly(zz,6,pz);
3096	lpoly(ell,9,pl);
3097
3098	for(i = 0; i < 4; i++) { //do 30 i= 1,4
3099	for(j = 0; j < 6; j++) { //do 30 j=1,6
3100	//elmax = elmax + emxcof[j][i]pz[j]pa[i];
3101	// elmax = elmax + emxcof[j][i]pz[i]pa[j];
3102	elmax = elmax + emxcof[i][j]pz[j]pa[i];
3103	}
3104	}
3105
3106	// assert(isnan(elmax) == false);
3107	(*selmax)=elmax;
3108
3109	// value of barrier at ang.mom. l
3110	if(il < 1){
3111	return;
3112	}
3113
3114	x = sel20/(*selmax);
3115	// assert(isnan(x) == false);
3116	y = sel80/(*selmax);
3117	// assert(isnan(y) == false);
3118
3119	if(el <= sel20) {
3120	// low l
3121	q = 0.2/(std::pow(sel20,2)std::pow(sel80,2)(sel20-sel80));
3122	qa = q(4.0std::pow(sel80,3) - std::pow(sel20,3));
3123	qb = -q(4.0std::pow(sel80,2) - std::pow(sel20,2));
3124	bfis = bfis(1.0 + qastd::pow(el,2) + qb*std::pow(el,3));
3125	}
3126	else {
3127	// high l
3128	aj = (-20.0std::pow(x,5) + 25.e0std::pow(x,4) - 4.0)std::pow((y-1.0),2)y*y;
3129	ak = (-20.0std::pow(y,5) + 25.0std::pow(y,4) - 1.0) * std::pow((x-1.0),2)xx;
3130	q = 0.2/(std::pow((y-x)((1.0-x)(1.0-y)xy),2));
3131	qa = q(ajy - ak*x);
3132	qb = -q(aj(2.0y + 1.0) - ak(2.0*x + 1.0));
3133	z = el/(*selmax);
3134	a1 = 4.0std::pow(z,5) - 5.0std::pow(z,4) + 1.0;
3135	a2 = qa(2.e0z + 1.e0);
3136	bfis=bfis(a1 + (z - 1.e0)(a2 + qbz)zz(z - 1.e0));
3137	}
3138
3139	if(bfis <= 0.0) {
3140	bfis=0.0;
3141	}
3142
3143	if(el > (*selmax)) {
3144	bfis=0.0;
3145	}
3146	(*sbfis)=bfis;
3147
3148	// now calculate rotating ground state energy
3149	if(el > (*selmax)) {
3150	return;
3151	}
3152
3153	for(k = 0; k < 4; k++) {
3154	for(l = 0; l < 6; l++) {
3155	for(m = 0; m < 5; m++) {
3156	//egs = egs + egscof[l][m][k]pz[l]pa[k]pl[2m-1];
3157	egs = egs + egscof[l][m][k]pz[l]pa[k]pl[2m];
3158	// egs = egs + egscof[m][l][k]pz[l]pa[k]pl[2m-1];
3159	}
3160	}
3161	}
3162
3163	(*segs)=egs;
3164	if((*segs) < 0.0) {
3165	(*segs)=0.0;
3166	}
3167
3168	return;
3169
3170	barfit900: //continue
3171	(*sbfis)=0.0;
3172	// for z<19 sbfis set to 1.0e3
3173	if (iz < 19) {
3174	(*sbfis) = 1.0e3;
3175	}
3176	(*segs)=0.0;
3177	(*selmax)=0.0;
3178	return;
3179
3180	barfit902:
3181	(*sbfis)=0.0;
3182	(*segs)=0.0;
3183	(*selmax)=0.0;
3184	return;
3185
3186	barfit910:
3187	(*sbfis)=0.0;
3188	(*segs)=0.0;
3189	(*selmax)=0.0;
3190	return;
3191
3192	barfit920:
3193	(*sbfis)=0.0;
3194	(*segs)=0.0;
3195	(*selmax)=0.0;
3196	return;
3197	}
3198
3199	G4double G4Abla::expohaz(G4int k, G4double T)
3200	{
3201	// TIRAGE ALEATOIRE DANS UNE EXPONENTIELLLE : Y=EXP(-X/T)
3202
3203	// assert(isnan((-1Tstd::log(haz(k)))) == false);
3204	return (-1.0Tstd::log(haz(k)));
3205	}
3206
3207	G4double G4Abla::fd(G4double E)
3208	{
3209	// DISTRIBUTION DE MAXWELL
3210
3211	return (E*std::exp(-E));
3212	}
3213
3214	G4double G4Abla::f(G4double E)
3215	{
3216	// FONCTION INTEGRALE DE FD(E)
3217	return (1.0 - (E + 1.0) * std::exp(-E));
3218	}
3219
3220	G4double G4Abla::fmaxhaz(G4double T)
3221	{
3222	// tirage aleatoire dans une maxwellienne
3223	// t : temperature
3224	//
3225	// declaration des variables
3226	//
3227
3228	const int pSize = 101;
3229	static G4double p[pSize];
3230
3231	// ial generateur pour le cascade (et les iy pour eviter les correlations)
3232	static G4int i = 0;
3233	static G4int itest = 0;
3234	// programme principal
3235
3236	// calcul des p(i) par approximation de newton
3237	p[pSize-1] = 8.0;
3238	G4double x = 0.1;
3239	G4double x1 = 0.0;
3240	G4double y = 0.0;
3241
3242	if (itest == 1) {
3243	goto fmaxhaz120;
3244	}
3245
3246	for(i = 1; i <= 99; i++) {
3247	fmaxhaz20:
3248	x1 = x - (f(x) - double(i)/100.0)/fd(x);
3249	x = x1;
3250	if (std::fabs(f(x) - double(i)/100.0) < 1e-5) {
3251	goto fmaxhaz100;
3252	}
3253	goto fmaxhaz20;
3254	fmaxhaz100:
3255	p[i] = x;
3256	} //end do
3257
3258	// itest = 1;
3259	itest = 0;
3260	// tirage aleatoire et calcul du x correspondant
3261	// par regression lineaire
3262	fmaxhaz120:
3263	standardRandom(&y, &(hazard->igraine[17]));
3264	i = nint(y*100);
3265
3266	// 2590 c ici on evite froidement les depassements de tableaux....(a.b. 3/9/99)
3267	if(i == 0) {
3268	goto fmaxhaz120;
3269	}
3270
3271	if (i == 1) {
3272	x = p[i]y100;
3273	}
3274	else {
3275	x = (p[i] - p[i-1])(y100 - i) + p[i];
3276	}
3277
3278	return(x*T);
3279	}
3280
3281	G4double G4Abla::pace2(G4double a, G4double z)
3282	{
3283	// PACE2
3284	// Cette fonction retourne le defaut de masse du noyau A,Z en MeV
3285	// RÃ©visÃ©e pour a, z flottants 25/4/2002 =
3286
3287	G4double pace2 = 0.0;
3288
3289	G4int ii = idint(a+0.5);
3290	G4int jj = idint(z+0.5);
3291
3292	if(ii <= 0 \|\| jj < 0) {
3293	pace2=0.;
3294	return pace2;
3295	}
3296
3297	if(jj > 300) {
3298	pace2=0.0;
3299	}
3300	else {
3301	pace2=pace->dm[ii][jj];
3302	}
3303	pace2=pace2/1000.;
3304
3305	if(pace->dm[ii][jj] == 0.) {
3306	if(ii < 12) {
3307	pace2=-500.;
3308	}
3309	else {
3310	guet(&a, &z, &pace2);
3311	pace2=pace2-ii*931.5;
3312	pace2=pace2/1000.;
3313	}
3314	}
3315
3316	return pace2;
3317	}
3318
3319	void G4Abla::guet(G4double x_par, G4double z_par, G4double *find_par)
3320	{
3321	// TABLE DE MASSES ET FORMULE DE MASSE TIRE DU PAPIER DE BRACK-GUET
3322	// Gives the theoritical value for mass excess...
3323	// RÃ©visÃ©e pour x, z flottants 25/4/2002
3324
3325	//real*8 x,z
3326	// dimension q(0:50,0:70)
3327	G4double x = (*x_par);
3328	G4double z = (*z_par);
3329	G4double find = (*find_par);
3330
3331	const G4int qrows = 50;
3332	const G4int qcols = 70;
3333	G4double q[qrows][qcols];
3334	for(G4int init_i = 0; init_i < qrows; init_i++) {
3335	for(G4int init_j = 0; init_j < qcols; init_j++) {
3336	q[init_i][init_j] = 0.0;
3337	}
3338	}
3339
3340	G4int ix=G4int(std::floor(x+0.5));
3341	G4int iz=G4int(std::floor(z+0.5));
3342	G4double zz = iz;
3343	G4double xx = ix;
3344	find = 0.0;
3345	G4double avol = 15.776;
3346	G4double asur = -17.22;
3347	G4double ac = -10.24;
3348	G4double azer = 8.0;
3349	G4double xjj = -30.03;
3350	G4double qq = -35.4;
3351	G4double c1 = -0.737;
3352	G4double c2 = 1.28;
3353
3354	if(ix <= 7) {
3355	q[0][1]=939.50;
3356	q[1][1]=938.21;
3357	q[1][2]=1876.1;
3358	q[1][3]=2809.39;
3359	q[2][4]=3728.34;
3360	q[2][3]=2809.4;
3361	q[2][5]=4668.8;
3362	q[2][6]=5606.5;
3363	q[3][5]=4669.1;
3364	q[3][6]=5602.9;
3365	q[3][7]=6535.27;
3366	q[4][6]=5607.3;
3367	q[4][7]=6536.1;
3368	q[5][7]=6548.3;
3369	find=q[iz][ix];
3370	}
3371	else {
3372	G4double xneu=xx-zz;
3373	G4double si=(xneu-zz)/xx;
3374	G4double x13=std::pow(xx,.333);
3375	G4double ee1=c1zzzz/x13;
3376	G4double ee2=c2zzzz/xx;
3377	G4double aux=1.+(9.*xjj/4./qq/x13);
3378	G4double ee3=xjjxxsi*si/aux;
3379	G4double ee4=avolxx+asur(std::pow(xx,.666))+ac*x13+azer;
3380	G4double tota = ee1 + ee2 + ee3 + ee4;
3381	find = 939.55xneu+938.77zz - tota;
3382	}
3383
3384	(*x_par) = x;
3385	(*z_par) = z;
3386	(*find_par) = find;
3387	}
3388
3389
3390	// Fission code
3391
3392	void G4Abla::even_odd(G4double r_origin,G4double r_even_odd,G4int &i_out)
3393	{
3394	// Procedure to calculate I_OUT from R_IN in a way that
3395	// on the average a flat distribution in R_IN results in a
3396	// fluctuating distribution in I_OUT with an even-odd effect as
3397	// given by R_EVEN_ODD
3398
3399	// /* ------------------------------------------------------------ */
3400	// /* EXAMPLES : */
3401	// /* ------------------------------------------------------------ */
3402	// /* If R_EVEN_ODD = 0 : */
3403	// /* CEIL(R_IN) ---- */
3404	// /* */
3405	// /* R_IN -> */
3406	// /* (somewhere in between CEIL(R_IN) and FLOOR(R_IN)) / /
3407	// /* */
3408	// /* FLOOR(R_IN) ---- --> I_OUT */
3409	// /* ------------------------------------------------------------ */
3410	// /* If R_EVEN_ODD > 0 : */
3411	// /* The interval for the above treatment is */
3412	// /* larger for FLOOR(R_IN) = even and */
3413	// /* smaller for FLOOR(R_IN) = odd */
3414	// /* For R_EVEN_ODD < 0 : just opposite treatment */
3415	// /* ------------------------------------------------------------ */
3416
3417	// /* ------------------------------------------------------------ */
3418	// /* On input: R_ORIGIN nuclear charge (real number) */
3419	// /* R_EVEN_ODD requested even-odd effect */
3420	// /* Intermediate quantity: R_IN = R_ORIGIN + 0.5 */
3421	// /* On output: I_OUT nuclear charge (integer) */
3422	// /* ------------------------------------------------------------ */
3423
3424	// G4double R_ORIGIN,R_IN,R_EVEN_ODD,R_REST,R_HELP;
3425	G4double r_in = 0.0, r_rest = 0.0, r_help = 0.0;
3426	G4double r_floor = 0.0;
3427	G4double r_middle = 0.0;
3428	// G4int I_OUT,N_FLOOR;
3429	G4int n_floor = 0;
3430
3431	r_in = r_origin + 0.5;
3432	r_floor = (float)((int)(r_in));
3433	if (r_even_odd < 0.001) {
3434	i_out = (int)(r_floor);
3435	}
3436	else {
3437	r_rest = r_in - r_floor;
3438	r_middle = r_floor + 0.5;
3439	n_floor = (int)(r_floor);
3440	if (n_floor%2 == 0) {
3441	// even before modif.
3442	r_help = r_middle + (r_rest - 0.5) * (1.0 - r_even_odd);
3443	}
3444	else {
3445	// odd before modification
3446	r_help = r_middle + (r_rest - 0.5) * (1.0 + r_even_odd);
3447	}
3448	i_out = (int)(r_help);
3449	}
3450	}
3451
3452	G4double G4Abla::umass(G4double z,G4double n,G4double beta)
3453	{
3454	// liquid-drop mass, Myers & Swiatecki, Lysekil, 1967
3455	// pure liquid drop, without pairing and shell effects
3456
3457	// On input: Z nuclear charge of nucleus
3458	// N number of neutrons in nucleus
3459	// beta deformation of nucleus
3460	// On output: binding energy of nucleus
3461
3462	G4double a = 0.0, umass = 0.0;
3463	G4double alpha = 0.0;
3464	G4double xcom = 0.0, xvs = 0.0, xe = 0.0;
3465	const G4double pi = 3.1416;
3466
3467	a = n + z;
3468	alpha = ( std::sqrt(5.0/(4.0pi)) ) beta;
3469	// assert(isnan(alpha) == false);
3470
3471	xcom = 1.0 - 1.7826 * ((a - 2.0z)/a)((a - 2.0*z)/a);
3472	// assert(isnan(xcom) == false);
3473	// factor for asymmetry dependence of surface and volume term
3474	xvs = - xcom * ( 15.4941 * a -
3475	17.9439 * std::pow(a,0.66667) * (1.0+0.4alphaalpha) );
3476	// sum of volume and surface energy
3477	xe = zz (0.7053/(std::pow(a,0.33333)) * (1.0-0.2alphaalpha) - 1.1529/a);
3478	// assert(isnan(xe) == false);
3479	umass = xvs + xe;
3480
3481	return umass;
3482	}
3483
3484	G4double G4Abla::ecoul(G4double z1,G4double n1,G4double beta1,G4double z2,G4double n2,G4double beta2,G4double d)
3485	{
3486	// Coulomb potential between two nuclei
3487	// surfaces are in a distance of d
3488	// in a tip to tip configuration
3489
3490	// approximate formulation
3491	// On input: Z1 nuclear charge of first nucleus
3492	// N1 number of neutrons in first nucleus
3493	// beta1 deformation of first nucleus
3494	// Z2 nuclear charge of second nucleus
3495	// N2 number of neutrons in second nucleus
3496	// beta2 deformation of second nucleus
3497	// d distance of surfaces of the nuclei
3498
3499	// G4double Z1,N1,beta1,Z2,N2,beta2,d,ecoul;
3500	G4double ecoul = 0;
3501	G4double dtot = 0;
3502	const G4double r0 = 1.16;
3503
3504	dtot = r0 * ( std::pow((z1+n1),0.33333) * (1.0+(2.0/3.0)*beta1)
3505	+ std::pow((z2+n2),0.33333) * (1.0+(2.0/3.0)*beta2) ) + d;
3506	ecoul = z1 * z2 * 1.44 / dtot;
3507
3508	// assert(isnan(ecoul) == false);
3509	return ecoul;
3510	}
3511
3512	void G4Abla::fissionDistri(G4double &a,G4double &z,G4double &e,
3513	G4double &a1,G4double &z1,G4double &e1,G4double &v1,
3514	G4double &a2,G4double &z2,G4double &e2,G4double &v2)
3515	{
3516	// On input: A, Z, E (mass, atomic number and exc. energy of compound nucleus
3517	// before fission)
3518	// On output: Ai, Zi, Ei (mass, atomic number and exc. energy of fragment 1 and 2
3519	// after fission)
3520
3521	// Additionally calculated but not put in the parameter list:
3522	// Kinetic energy of prefragments EkinR1, EkinR2
3523
3524	// Translation of SIMFIS18.PLI (KHS, 2.1.2001)
3525
3526	// This program calculates isotopic distributions of fission fragments
3527	// with a semiempirical model
3528	// Copy from SIMFIS3, KHS, 8. February 1995
3529	// Modifications made by Jose Benlliure and KHS in August 1996
3530	// Energy counted from lowest barrier (J. Benlliure, KHS 1997)
3531	// Some bugs corrected (J. Benlliure, KHS 1997)
3532	// Version used for thesis S. Steinhaueser (August 1997)
3533	// (Curvature of LD potential increased by factor of 2!)
3534
3535	// Weiter veraendert mit der Absicht, eine Version zu erhalten, die
3536	// derjenigen entspricht, die von J. Benlliure et al.
3537	// in Nucl. Phys. A 628 (1998) 458 verwendet wurde,
3538	// allerdings ohne volle Neutronenabdampfung.
3539
3540	// The excitation energy was calculate now for each fission channel
3541	// separately. The dissipation from saddle to scission was taken from
3542	// systematics, the deformation energy at scission considers the shell
3543	// effects in a simplified way, and the fluctuation is included.
3544	// KHS, April 1999
3545
3546	// The width in N/Z was carefully adapted to values given by Lang et al.
3547
3548	// The width and eventually a shift in N/Z (polarization) follows the
3549	// following rules:
3550
3551	// The line N/Z following UCD has an angle of std::atan(Zcn/Ncn)
3552	// to the horizontal axis on a chart of nuclides.
3553	// (For 238U the angle is 32.2 deg.)
3554
3555	// The following relations hold: (from Armbruster)
3556	//
3557	// sigma(N) (A=const) = sigma(Z) (A=const)
3558	// sigma(A) (N=const) = sigma(Z) (N=const)
3559	// sigma(A) (Z=const) = sigma(N) (Z=const)
3560	//
3561	// From this we get:
3562	// sigma(Z) (N=const) * N = sigma(N) (Z=const) * Z
3563	// sigma(A) (Z=const) = sigma(Z) (A=const) * A/Z
3564	// sigma(N) (Z=const) = sigma(Z) (A=const) * A/Z
3565	// Zsigma(N) (Z=const) = Nsigma(Z) (N=const) = A*sigma(Z) (A=const)
3566
3567	// Excitation energy now calculated above the lowest potential point
3568	// Inclusion of a distribution of excitation energies
3569
3570	// Several modifications, starting from SIMFIS12: KHS November 2000
3571	// This version seems to work quite well for 238U.
3572	// The transition from symmetric to asymmetric fission around 226Th
3573	// is reasonably well reproduced, although St. I is too strong and St. II
3574	// is too weak. St. I and St. II are also weakly seen for 208Pb.
3575
3576	// Extensions for an event generator of fission events (21.11.2000,KHS)
3577
3578	// Defalt parameters (IPARS) rather carefully adjusted to
3579	// pre-neutron mass distributions of Vives et al. (238U + n)
3580	// Die Parameter Fgamma1 und Fgamma2 sind kleiner als die resultierenden
3581	// Breiten der Massenverteilungen!!!
3582	// Fgamma1 und Fgamma2 wurden angepaï¿œ, so daï¿œ
3583	// Sigma-A(ST-I) = 3.3, Sigma-A(St-II) = 5.8 (nach Vives)
3584
3585	// Parameters of the model carefully adjusted by KHS (2.2.2001) to
3586	// 238U + 208Pb, 1000 A MeV, Timo Enqvist et al.
3587
3588
3589	G4double n = 0.0;
3590	G4double nlight1 = 0.0, nlight2 = 0.0;
3591	G4double aheavy1 = 0.0,alight1 = 0.0, aheavy2 = 0.0, alight2 = 0.0;
3592	G4double eheavy1 = 0.0, elight1 = 0.0, eheavy2 = 0.0, elight2 = 0.0;
3593	G4double zheavy1_shell = 0.0, zheavy2_shell = 0.0;
3594	G4double zlight1 = 0.0, zlight2 = 0.0;
3595	G4double masscurv = 0.0;
3596	G4double sasymm1 = 0.0, sasymm2 = 0.0, ssymm = 0.0, ysum = 0.0, yasymm = 0.0;
3597	G4double ssymm_mode1 = 0.0, ssymm_mode2 = 0.0;
3598	G4double cz_asymm1_saddle = 0.0, cz_asymm2_saddle = 0.0;
3599	// Curvature at saddle, modified by ld-potential
3600	G4double wzasymm1_saddle, wzasymm2_saddle, wzsymm_saddle = 0.0;
3601	G4double wzasymm1_scission = 0.0, wzasymm2_scission = 0.0, wzsymm_scission = 0.0;
3602	G4double wzasymm1 = 0.0, wzasymm2 = 0.0, wzsymm = 0.0;
3603	G4double nlight1_eff = 0.0, nlight2_eff = 0.0;
3604	G4int imode = 0;
3605	G4double rmode = 0.0;
3606	G4double z1mean = 0.0, z2mean = 0.0, z1width = 0.0, za1width = 0.0;
3607	// G4double Z1,Z2,N1R,N2R,A1R,A2R,N1,N2,A1,A2;
3608	G4double n1r = 0.0, n2r = 0.0, a1r = 0.0, a2r = 0.0, n1 = 0.0, n2 = 0.0;
3609
3610	G4double zsymm = 0.0, nsymm = 0.0, asymm = 0.0;
3611	G4double n1mean = 0.0, n2mean, n1width;
3612	G4double dueff = 0.0;
3613	// effective shell effect at lowest barrier
3614	G4double eld = 0.0;
3615	// Excitation energy with respect to ld barrier
3616	G4double re1 = 0.0, re2 = 0.0, re3 = 0.0;
3617	G4double eps1 = 0.0, eps2 = 0.0;
3618	G4double n1ucd = 0.0, n2ucd = 0.0, z1ucd = 0.0, z2ucd = 0.0;
3619	G4double beta = 0.0, beta1 = 0.0, beta2 = 0.0;
3620
3621	G4double dn1_pol = 0.0;
3622	// shift of most probable neutron number for given Z,
3623	// according to polarization
3624	G4int i_help = 0;
3625
3626	// /* Parameters of the semiempirical fission model */
3627	G4double a_levdens = 0.0;
3628	// /* level-density parameter */
3629	G4double a_levdens_light1 = 0.0, a_levdens_light2 = 0.0;
3630	G4double a_levdens_heavy1 = 0.0, a_levdens_heavy2 = 0.0;
3631	const G4double r_null = 1.16;
3632	// /* radius parameter */
3633	G4double epsilon_1_saddle = 0.0, epsilon0_1_saddle = 0.0;
3634	G4double epsilon_2_saddle = 0.0, epsilon0_2_saddle = 0.0, epsilon_symm_saddle = 0.0;
3635	G4double epsilon_1_scission = 0.0, epsilon0_1_scission = 0.0;
3636	G4double epsilon_2_scission = 0.0, epsilon0_2_scission = 0.0;
3637	G4double epsilon_symm_scission = 0.0;
3638	// /* modified energy */
3639	G4double e_eff1_saddle = 0.0, e_eff2_saddle = 0.0;
3640	G4double epot0_mode1_saddle = 0.0, epot0_mode2_saddle = 0.0, epot0_symm_saddle = 0.0;
3641	G4double epot_mode1_saddle = 0.0, epot_mode2_saddle = 0.0, epot_symm_saddle = 0.0;
3642	G4double e_defo = 0.0, e_defo1 = 0.0, e_defo2 = 0.0, e_scission = 0.0, e_asym = 0.0;
3643	G4double e1exc = 0.0, e2exc = 0.0;
3644	G4double e1exc_sigma = 0.0, e2exc_sigma = 0.0;
3645	G4double e1final = 0.0, e2final = 0.0;
3646
3647	const G4double r0 = 1.16;
3648	G4double tker = 0.0;
3649	G4double ekin1 = 0.0, ekin2 = 0.0;
3650	// G4double EkinR1,EkinR2,E1,E2,V1,V2;
3651	G4double ekinr1 = 0.0, ekinr2 = 0.0;
3652	G4int icz = 0, k = 0;
3653
3654	// Input parameters:
3655	//OMMENT(Nuclear charge number);
3656	// G4double Z;
3657	//OMMENT(Nuclear mass number);
3658	// G4double A;
3659	//OMMENT(Excitation energy above fission barrier);
3660	// G4double E;
3661
3662	// Model parameters:
3663	//OMMENT(position of heavy peak valley 1);
3664	const G4double nheavy1 = 83.0;
3665	//OMMENT(position of heavy peak valley 2);
3666	const G4double nheavy2 = 90.0;
3667	//OMMENT(Shell effect for valley 1);
3668	const G4double delta_u1_shell = -2.65;
3669	// Parameter (Delta_U1_shell = -2)
3670	//OMMENT(Shell effect for valley 2);
3671	const G4double delta_u2_shell = -3.8;
3672	// Parameter (Delta_U2_shell = -3.2)
3673	//OMMENT(I: used shell effect);
3674	G4double delta_u1 = 0.0;
3675	//omment(I: used shell effect);
3676	G4double delta_u2 = 0.0;
3677	//OMMENT(Curvature of asymmetric valley 1);
3678	const G4double cz_asymm1_shell = 0.7;
3679	//OMMENT(Curvature of asymmetric valley 2);
3680	const G4double cz_asymm2_shell = 0.15;
3681	//OMMENT(Factor for width of distr. valley 1);
3682	const G4double fwidth_asymm1 = 0.63;
3683	//OMMENT(Factor for width of distr. valley 2);
3684	const G4double fwidth_asymm2 = 0.97;
3685	// Parameter (CZ_asymm2_scission = 0.12)
3686	//OMMENT(Parameter x: a = A/x);
3687	const G4double xlevdens = 12.0;
3688	//OMMENT(Factor to gamma_heavy1);
3689	const G4double fgamma1 = 2.0;
3690	//OMMENT(I: fading of shells (general));
3691	G4double gamma = 0.0;
3692	//OMMENT(I: fading of shell 1);
3693	G4double gamma_heavy1 = 0.0;
3694	//OMMENT(I: fading of shell 2);
3695	G4double gamma_heavy2 = 0.0;
3696	//OMMENT(Zero-point energy at saddle);
3697	const G4double e_zero_point = 0.5;
3698	//OMMENT(I: friction from saddle to scission);
3699	G4double e_saddle_scission = 0.0;
3700	//OMMENT(Friction factor);
3701	const G4double friction_factor = 1.0;
3702	//OMMENT(I: Internal counter for different modes); INIT(0,0,0)
3703	// Integer*4 I_MODE(3)
3704	//OMMENT(I: Yield of symmetric mode);
3705	G4double ysymm = 0.0;
3706	//OMMENT(I: Yield of asymmetric mode 1);
3707	G4double yasymm1 = 0.0;
3708	//OMMENT(I: Yield of asymmetric mode 2);
3709	G4double yasymm2 = 0.0;
3710	//OMMENT(I: Effective position of valley 1);
3711	G4double nheavy1_eff = 0.0;
3712	//OMMENT(I: position of heavy peak valley 1);
3713	G4double zheavy1 = 0.0;
3714	//omment(I: Effective position of valley 2);
3715	G4double nheavy2_eff = 0.0;
3716	//OMMENT(I: position of heavy peak valley 2);
3717	G4double zheavy2 = 0.0;
3718	//omment(I: Excitation energy above saddle 1);
3719	G4double eexc1_saddle = 0.0;
3720	//omment(I: Excitation energy above saddle 2);
3721	G4double eexc2_saddle = 0.0;
3722	//omment(I: Excitation energy above lowest saddle);
3723	G4double eexc_max = 0.0;
3724	//omment(I: Effective mass mode 1);
3725	G4double aheavy1_mean = 0.0;
3726	//omment(I: Effective mass mode 2);
3727	G4double aheavy2_mean = 0.0;
3728	//omment(I: Width of symmetric mode);
3729	G4double wasymm_saddle = 0.0;
3730	//OMMENT(I: Width of asymmetric mode 1);
3731	G4double waheavy1_saddle = 0.0;
3732	//OMMENT(I: Width of asymmetric mode 2);
3733	G4double waheavy2_saddle = 0.0;
3734	//omment(I: Width of symmetric mode);
3735	G4double wasymm = 0.0;
3736	//OMMENT(I: Width of asymmetric mode 1);
3737	G4double waheavy1 = 0.0;
3738	//OMMENT(I: Width of asymmetric mode 2);
3739	G4double waheavy2 = 0.0;
3740	//OMMENT(I: Even-odd effect in Z);
3741	G4double r_e_o = 0.0, r_e_o_exp = 0.0;
3742	//OMMENT(I: Curveture of symmetric valley);
3743	G4double cz_symm = 0.0;
3744	//OMMENT(I: Curvature of mass distribution for fixed Z);
3745	G4double cn = 0.0;
3746	//OMMENT(I: Curvature of Z distribution for fixed A);
3747	G4double cz = 0.0;
3748	//OMMENT(Minimum neutron width for constant Z);
3749	const G4double sigzmin = 1.16;
3750	//OMMENT(Surface distance of scission configuration);
3751	const G4double d = 2.0;
3752
3753	// /* Charge polarisation from Wagemanns p. 397: */
3754	//OMMENT(Charge polarisation standard I);
3755	const G4double cpol1 = 0.65;
3756	//OMMENT(Charge polarisation standard II);
3757	const G4double cpol2 = 0.55;
3758	//OMMENT(=1: Polarisation simult. in N and Z);
3759	const G4int nzpol = 1;
3760	//OMMENT(=1: test output, =0: no test output);
3761	const G4int itest = 0;
3762
3763	// G4double UMASS, ECOUL, reps1, reps2, rn1_pol;
3764	G4double reps1 = 0.0, reps2 = 0.0, rn1_pol = 0.0;
3765	// Float_t HAZ,GAUSSHAZ;
3766	G4int kkk = 0;
3767	// G4int kkk = 10; // PK
3768
3769	// I_MODE = 0;
3770
3771	if(itest == 1) {
3772	G4cout << " cn mass " << a << G4endl;
3773	G4cout << " cn charge " << z << G4endl;
3774	G4cout << " cn energy " << e << G4endl;
3775	}
3776
3777	// /* average Z of asymmetric and symmetric components: */
3778	n = a - z; /* neutron number of the fissioning nucleus */
3779
3780	k = 0;
3781	icz = 0;
3782	if ( (std::pow(z,2)/a < 25.0) \|\| (n < nheavy2) \|\| (e > 500.0) ) {
3783	icz = -1;
3784	// GOTO 1002;
3785	goto milledeux;
3786	}
3787
3788	nlight1 = n - nheavy1;
3789	nlight2 = n - nheavy2;
3790
3791	// /* Polarisation assumed for standard I and standard II:
3792	// Z - Zucd = cpol (for A = const);
3793	// from this we get (see Armbruster)
3794	// Z - Zucd = Acn/Ncn * cpol (for N = const) */
3795
3796	zheavy1_shell = ((nheavy1/n) * z) - ((a/n) * cpol1);
3797	zheavy2_shell = ((nheavy2/n) * z) - ((a/n) * cpol2);
3798
3799	e_saddle_scission =
3800	(-24.0 + 0.02227 * (std::pow(z,2))/(std::pow(a,0.33333)) ) * friction_factor;
3801
3802	// /* Energy dissipated from saddle to scission */
3803	// /* F. Rejmund et al., Nucl. Phys. A 678 (2000) 215, fig. 4 b */
3804	// E_saddle_scission = DMAX1(0.,E_saddle_scission);
3805	if (e_saddle_scission > 0.) {
3806	e_saddle_scission = e_saddle_scission;
3807	}
3808	else {
3809	e_saddle_scission = 0.;
3810	}
3811	// /* Semiempirical fission model: */
3812
3813	// /* Fit to experimental result on curvature of potential at saddle */
3814	// /* reference: */
3815	// /* IF Z**2/A < 33.15E0 THEN
3816	// MassCurv = 30.5438538E0 - 4.00212049E0 * Z**2/A
3817	// + 0.11983384E0 * Z4 / (A2) ;
3818	// ELSE
3819	// MassCurv = 10.E0 ** (7.16993332E0 - 0.26602401E0 * Z**2/A
3820	// + 0.00283802E0 * Z4 / (A2)) ; */
3821	// /* New parametrization of T. Enqvist according to Mulgin et al. 1998 */
3822	if ( (std::pow(z,2))/a < 34.0) {
3823	masscurv = std::pow( 10.0,(-1.093364 + 0.082933 * (std::pow(z,2)/a)
3824	- 0.0002602 * (std::pow(z,4)/std::pow(a,2))) );
3825	} else {
3826	masscurv = std::pow( 10.0,(3.053536 - 0.056477 * (std::pow(z,2)/a)
3827	+ 0.0002454 * (std::pow(z,4)/std::pow(a,2))) );
3828	}
3829
3830	cz_symm = (8.0/std::pow(z,2)) * masscurv;
3831
3832	if(itest == 1) {
3833	G4cout << "cz_symmetry= " << cz_symm << G4endl;
3834	}
3835
3836	if (cz_symm < 0) {
3837	icz = -1;
3838	// GOTO 1002;
3839	goto milledeux;
3840	}
3841
3842	// /* proton number in symmetric fission (centre) */
3843	zsymm = z/2.0;
3844	nsymm = n/2.0;
3845	asymm = nsymm + zsymm;
3846
3847	zheavy1 = (cz_symmzsymm + cz_asymm1_shellzheavy1_shell)/(cz_symm + cz_asymm1_shell);
3848	zheavy2 = (cz_symmzsymm + cz_asymm2_shellzheavy2_shell)/(cz_symm + cz_asymm2_shell);
3849	// /* position of valley due to influence of liquid-drop potential */
3850	nheavy1_eff = (zheavy1 + (a/n * cpol1))*(n/z);
3851	nheavy2_eff = (zheavy2 + (a/n * cpol2))*(n/z);
3852	nlight1_eff = n - nheavy1_eff;
3853	nlight2_eff = n - nheavy2_eff;
3854	// /* proton number of light fragments (centre) */
3855	zlight1 = z - zheavy1;
3856	// /* proton number of light fragments (centre) */
3857	zlight2 = z - zheavy2;
3858	aheavy1 = nheavy1_eff + zheavy1;
3859	aheavy2 = nheavy2_eff + zheavy2;
3860	aheavy1_mean = aheavy1;
3861	aheavy2_mean = aheavy2;
3862	alight1 = nlight1_eff + zlight1;
3863	alight2 = nlight2_eff + zlight2;
3864
3865	a_levdens = a / xlevdens;
3866	a_levdens_heavy1 = aheavy1 / xlevdens;
3867	a_levdens_heavy2 = aheavy2 / xlevdens;
3868	a_levdens_light1 = alight1 / xlevdens;
3869	a_levdens_light2 = alight2 / xlevdens;
3870	gamma = a_levdens / (0.4 * (std::pow(a,1.3333)) );
3871	gamma_heavy1 = ( a_levdens_heavy1 / (0.4 * (std::pow(aheavy1,1.3333)) ) ) * fgamma1;
3872	gamma_heavy2 = a_levdens_heavy2 / (0.4 * (std::pow(aheavy2,1.3333)) );
3873
3874	cz_asymm1_saddle = cz_asymm1_shell + cz_symm;
3875	cz_asymm2_saddle = cz_asymm2_shell + cz_symm;
3876
3877	// Up to here: Ok! Checked CS 10/10/05
3878
3879	cn = umass(zsymm,(nsymm+1.),0.0) + umass(zsymm,(nsymm-1.),0.0)
3880	+ 1.44 * (std::pow(zsymm,2))/
3881	( (std::pow(r_null,2)) *
3882	( std::pow((asymm+1.0),0.33333) + std::pow((asymm-1.0),0.33333) ) *
3883	( std::pow((asymm+1.0),0.33333) + std::pow((asymm-1.0),0.33333) ) )
3884	- 2.0 * umass(zsymm,nsymm,0.0)
3885	- 1.44 * (std::pow(zsymm,2))/
3886	( ( 2.0 * r_null * (std::pow(asymm,0.33333)) ) *
3887	( 2.0 * r_null * (std::pow(asymm,0.33333)) ) );
3888
3889	// /* shell effect in valley of mode 1 */
3890	delta_u1 = delta_u1_shell + (std::pow((zheavy1_shell-zheavy1),2))*cz_asymm1_shell;
3891	// /* shell effect in valley of mode 2 */
3892	delta_u2 = delta_u2_shell + (std::pow((zheavy2_shell-zheavy2),2))*cz_asymm2_shell;
3893
3894	// /* liquid drop energies
3895	// at the centres of the different shell effects
3896	// with respect to liquid drop at symmetry: */
3897	epot0_mode1_saddle = (std::pow((zheavy1-zsymm),2)) * cz_symm;
3898	epot0_mode2_saddle = (std::pow((zheavy2-zsymm),2)) * cz_symm;
3899	epot0_symm_saddle = 0.0;
3900
3901	if (itest == 1) {
3902	G4cout << "check zheavy1 = " << zheavy1 << G4endl;
3903	G4cout << "check zheavy2 = " << zheavy2 << G4endl;
3904	G4cout << "check zsymm = " << zsymm << G4endl;
3905	G4cout << "check czsymm = " << cz_symm << G4endl;
3906	G4cout << "check epot0_mode1_saddle = " << epot0_mode1_saddle << G4endl;
3907	G4cout << "check epot0_mode2_saddle = " << epot0_mode2_saddle << G4endl;
3908	G4cout << "check epot0_symm_saddle = " << epot0_symm_saddle << G4endl;
3909	G4cout << "delta_u1 = " << delta_u1 << G4endl;
3910	G4cout << "delta_u2 = " << delta_u2 << G4endl;
3911	}
3912
3913	// /* energies including shell effects
3914	// at the centres of the different shell effects
3915	// with respect to liquid drop at symmetry: */
3916	epot_mode1_saddle = epot0_mode1_saddle + delta_u1;
3917	epot_mode2_saddle = epot0_mode2_saddle + delta_u2;
3918	epot_symm_saddle = epot0_symm_saddle;
3919	if (itest == 1) {
3920	G4cout << "check epot_mode1_saddle = " << epot_mode1_saddle << G4endl;
3921	G4cout << "check epot_mode2_saddle = " << epot_mode2_saddle << G4endl;
3922	G4cout << "check epot_symm_saddle = " << epot_symm_saddle << G4endl;
3923	}
3924
3925	// /* Minimum of potential with respect to ld potential at symmetry */
3926	dueff = min(epot_mode1_saddle,epot_mode2_saddle);
3927	dueff = min(dueff,epot_symm_saddle);
3928	dueff = dueff - epot_symm_saddle;
3929
3930	eld = e + dueff + e_zero_point;
3931
3932	if (itest == 1) {
3933	G4cout << "check dueff = " << dueff << G4endl;
3934	G4cout << "check e = " << e << G4endl;
3935	G4cout << "check e_zero_point = " << e_zero_point << G4endl;
3936	G4cout << "check eld = " << eld << G4endl;
3937	}
3938	// Up to here: Ok! Checked CS 10/10/05
3939
3940	// /* E = energy above lowest effective barrier */
3941	// /* Eld = energy above liquid-drop barrier */
3942
3943	// /* Due to this treatment the energy E on input means the excitation */
3944	// /* energy above the lowest saddle. */
3945
3946	// /* These energies are not used */
3947	eheavy1 = e * aheavy1 / a;
3948	eheavy2 = e * aheavy2 / a;
3949	elight1 = e * alight1 / a;
3950	elight2 = e * alight2 / a;
3951
3952	epsilon0_1_saddle = eld - e_zero_point - epot0_mode1_saddle;
3953	// /* excitation energy at saddle mode 1 without shell effect */
3954	epsilon0_2_saddle = eld - e_zero_point - epot0_mode2_saddle;
3955	// /* excitation energy at saddle mode 2 without shell effect */
3956
3957	epsilon_1_saddle = eld - e_zero_point - epot_mode1_saddle;
3958	// /* excitation energy at saddle mode 1 with shell effect */
3959	epsilon_2_saddle = eld - e_zero_point - epot_mode2_saddle;
3960	// /* excitation energy at saddle mode 2 with shell effect */
3961	epsilon_symm_saddle = eld - e_zero_point - epot_symm_saddle;
3962
3963	// /* global parameters */
3964	eexc1_saddle = epsilon_1_saddle;
3965	eexc2_saddle = epsilon_2_saddle;
3966	eexc_max = max(eexc1_saddle,eexc2_saddle);
3967	eexc_max = max(eexc_max,eld);
3968
3969	// /* EEXC_MAX is energy above the lowest saddle */
3970
3971
3972	epsilon0_1_scission = eld + e_saddle_scission - epot0_mode1_saddle;
3973	// /* excitation energy without shell effect */
3974	epsilon0_2_scission = eld + e_saddle_scission - epot0_mode2_saddle;
3975	// /* excitation energy without shell effect */
3976
3977	epsilon_1_scission = eld + e_saddle_scission - epot_mode1_saddle;
3978	// /* excitation energy at scission */
3979	epsilon_2_scission = eld+ e_saddle_scission - epot_mode2_saddle;
3980	// /* excitation energy at scission */
3981	epsilon_symm_scission = eld + e_saddle_scission - epot_symm_saddle;
3982	// /* excitation energy of symmetric fragment at scission */
3983
3984	// /* Calculate widhts at the saddle: */
3985
3986	e_eff1_saddle = epsilon0_1_saddle - delta_u1 * (std::exp((-epsilon_1_saddle*gamma)));
3987
3988	if (e_eff1_saddle > 0.0) {
3989	wzasymm1_saddle = std::sqrt( (0.5 *
3990	(std::sqrt(1.0/a_levdens*e_eff1_saddle)) /
3991	(cz_asymm1_shell * std::exp((-epsilon_1_saddle*gamma)) + cz_symm) ) );
3992	}
3993	else {
3994	wzasymm1_saddle = 1.0;
3995	}
3996
3997	e_eff2_saddle = epsilon0_2_saddle - delta_u2 * (std::exp((-epsilon_2_saddle*gamma)));
3998	if (e_eff2_saddle > 0.0) {
3999	wzasymm2_saddle = std::sqrt( (0.5 *
4000	(std::sqrt(1.0/a_levdens*e_eff2_saddle)) /
4001	(cz_asymm2_shell * std::exp((-epsilon_2_saddle*gamma)) + cz_symm) ) );
4002	}
4003	else {
4004	wzasymm2_saddle = 1.0;
4005	}
4006
4007	if (eld > e_zero_point) {
4008	if ( (eld + epsilon_symm_saddle) < 0.0) {
4009	G4cout << "<e> eld + epsilon_symm_saddle < 0" << G4endl;
4010	}
4011	wzsymm_saddle = std::sqrt( (0.5 *
4012	(std::sqrt(1.0/a_levdens*(eld+epsilon_symm_saddle))) / cz_symm ) );
4013	} else {
4014	wzsymm_saddle = 1.0;
4015	}
4016
4017	if (itest == 1) {
4018	G4cout << "wz1(saddle) = " << wzasymm1_saddle << G4endl;
4019	G4cout << "wz2(saddle) = " << wzasymm2_saddle << G4endl;
4020	G4cout << "wzsymm(saddle) = " << wzsymm_saddle << G4endl;
4021	}
4022
4023	// /* Calculate widhts at the scission point: */
4024	// /* fits of ref. Beizin 1991 (Plots brought to GSI by Sergei Zhdanov) */
4025
4026	wzsymm_scission = wzsymm_saddle;
4027
4028	if (e_saddle_scission == 0.0) {
4029
4030	wzasymm1_scission = wzasymm1_saddle;
4031	wzasymm2_scission = wzasymm2_saddle;
4032
4033	}
4034	else {
4035
4036	if (nheavy1_eff > 75.0) {
4037	wzasymm1_scission = (std::sqrt(21.0)) * z/a;
4038	wzasymm2_scission = (std::sqrt (max( (70.0-28.0)/3.0(zz/a-35.0)+28.,0.0 )) ) * z/a;
4039	}
4040	else {
4041	wzasymm1_scission = wzasymm1_saddle;
4042	wzasymm2_scission = wzasymm2_saddle;
4043	}
4044
4045	}
4046
4047	wzasymm1_scission = max(wzasymm1_scission,wzasymm1_saddle);
4048	wzasymm2_scission = max(wzasymm2_scission,wzasymm2_saddle);
4049
4050	wzasymm1 = wzasymm1_scission * fwidth_asymm1;
4051	wzasymm2 = wzasymm2_scission * fwidth_asymm2;
4052	wzsymm = wzsymm_scission;
4053
4054	/* if (ITEST == 1) {
4055	G4cout << "WZ1(scission) = " << WZasymm1_scission << G4endl;
4056	G4cout << "WZ2(scission) = " << WZasymm2_scission << G4endl;
4057	G4cout << "WZsymm(scission) = " << WZsymm_scission << G4endl;
4058	}
4059	if (ITEST == 1) {
4060	G4cout << "WZ1(scission) final= " << WZasymm1 << G4endl;
4061	G4cout << "WZ2(scission) final= " << WZasymm2 << G4endl;
4062	G4cout << "WZsymm(scission) final= " << WZsymm << G4endl;
4063	} */
4064
4065	wasymm = wzsymm * a/z;
4066	waheavy1 = wzasymm1 * a/z;
4067	waheavy2 = wzasymm2 * a/z;
4068
4069	wasymm_saddle = wzsymm_saddle * a/z;
4070	waheavy1_saddle = wzasymm1_saddle * a/z;
4071	waheavy2_saddle = wzasymm2_saddle * a/z;
4072
4073	if (itest == 1) {
4074	G4cout << "wasymm = " << wzsymm << G4endl;
4075	G4cout << "waheavy1 = " << waheavy1 << G4endl;
4076	G4cout << "waheavy2 = " << waheavy2 << G4endl;
4077	}
4078	// Up to here: Ok! Checked CS 11/10/05
4079
4080	if ( (epsilon0_1_saddle - delta_u1std::exp((-epsilon_1_saddlegamma_heavy1))) < 0.0) {
4081	sasymm1 = -10.0;
4082	}
4083	else {
4084	sasymm1 = 2.0 * std::sqrt( a_levdens * (epsilon0_1_saddle -
4085	delta_u1(std::exp((-epsilon_1_saddlegamma_heavy1))) ) );
4086	}
4087
4088	if ( (epsilon0_2_saddle - delta_u2std::exp((-epsilon_2_saddlegamma_heavy2))) < 0.0) {
4089	sasymm2 = -10.0;
4090	}
4091	else {
4092	sasymm2 = 2.0 * std::sqrt( a_levdens * (epsilon0_2_saddle -
4093	delta_u2(std::exp((-epsilon_2_saddlegamma_heavy2))) ) );
4094	}
4095
4096	if (epsilon_symm_saddle > 0.0) {
4097	ssymm = 2.0 * std::sqrt( a_levdens*(epsilon_symm_saddle) );
4098	}
4099	else {
4100	ssymm = -10.0;
4101	}
4102
4103	if (ssymm > -10.0) {
4104	ysymm = 1.0;
4105
4106	if (epsilon0_1_saddle < 0.0) {
4107	// /* low energy */
4108	yasymm1 = std::exp((sasymm1-ssymm)) * wzasymm1_saddle / wzsymm_saddle * 2.0;
4109	// /* factor of 2 for symmetry classes */
4110	}
4111	else {
4112	// /* high energy */
4113	ssymm_mode1 = 2.0 * std::sqrt( a_levdens*(epsilon0_1_saddle) );
4114	yasymm1 = ( std::exp((sasymm1-ssymm)) - std::exp((ssymm_mode1 - ssymm)) )
4115	* wzasymm1_saddle / wzsymm_saddle * 2.0;
4116	}
4117
4118	if (epsilon0_2_saddle < 0.0) {
4119	// /* low energy */
4120	yasymm2 = std::exp((sasymm2-ssymm)) * wzasymm2_saddle / wzsymm_saddle * 2.0;
4121	// /* factor of 2 for symmetry classes */
4122	}
4123	else {
4124	// /* high energy */
4125	ssymm_mode2 = 2.0 * std::sqrt( a_levdens*(epsilon0_2_saddle) );
4126	yasymm2 = ( std::exp((sasymm2-ssymm)) - std::exp((ssymm_mode2 - ssymm)) )
4127	* wzasymm2_saddle / wzsymm_saddle * 2.0;
4128	}
4129	// /* difference in the exponent in order */
4130	// /* to avoid numerical overflow */
4131
4132	}
4133	else {
4134	if ( (sasymm1 > -10.0) && (sasymm2 > -10.0) ) {
4135	ysymm = 0.0;
4136	yasymm1 = std::exp(sasymm1) * wzasymm1_saddle * 2.0;
4137	yasymm2 = std::exp(sasymm2) * wzasymm2_saddle * 2.0;
4138	}
4139	}
4140
4141	// /* normalize */
4142	ysum = ysymm + yasymm1 + yasymm2;
4143	if (ysum > 0.0) {
4144	ysymm = ysymm / ysum;
4145	yasymm1 = yasymm1 / ysum;
4146	yasymm2 = yasymm2 / ysum;
4147	yasymm = yasymm1 + yasymm2;
4148	}
4149	else {
4150	ysymm = 0.0;
4151	yasymm1 = 0.0;
4152	yasymm2 = 0.0;
4153	// /* search minimum threshold and attribute all events to this mode */
4154	if ( (epsilon_symm_saddle < epsilon_1_saddle) && (epsilon_symm_saddle < epsilon_2_saddle) ) {
4155	ysymm = 1.0;
4156	}
4157	else {
4158	if (epsilon_1_saddle < epsilon_2_saddle) {
4159	yasymm1 = 1.0;
4160	}
4161	else {
4162	yasymm2 = 1.0;
4163	}
4164	}
4165	}
4166
4167	if (itest == 1) {
4168	G4cout << "ysymm normalized= " << ysymm << G4endl;
4169	G4cout << "yasymm1 normalized= " << yasymm1 << G4endl;
4170	G4cout << "yasymm2 normalized= " << yasymm2 << G4endl;
4171	}
4172	// Up to here: Ok! Ckecked CS 11/10/05
4173
4174	// /* even-odd effect */
4175	// /* simple parametrization KHS, Nov. 2000. From Rejmund et al. */
4176	if ((int)(z) % 2 == 0) {
4177	r_e_o_exp = -0.017 * (e_saddle_scission + eld) * (e_saddle_scission + eld);
4178	if ( r_e_o_exp < -307.0) {
4179	r_e_o_exp = -307.0;
4180	r_e_o = std::pow(10.0,r_e_o_exp);
4181	}
4182	else {
4183	r_e_o = std::pow(10.0,r_e_o_exp);
4184	}
4185	}
4186	else {
4187	r_e_o = 0.0;
4188	}
4189
4190	// $LOOP; /* event loop */
4191	// I_COUNT = I_COUNT + 1;
4192
4193	// /* random decision: symmetric or asymmetric */
4194	// /* IMODE = 1 means asymmetric fission, mode 1,
4195	// IMODE = 2 means asymmetric fission, mode 2,
4196	// IMODE = 3 means symmetric */
4197	// RMODE = dble(HAZ(k));
4198	// rmode = rnd.rndm();
4199
4200	// Safety check added to make sure we always select well defined
4201	// fission mode.
4202	do {
4203	rmode = haz(k);
4204	// Cast for test CS 11/10/05
4205	// RMODE = 0.54;
4206	// rmode = 0.54;
4207	if (rmode < yasymm1) {
4208	imode = 1;
4209	} else if ( (rmode > yasymm1) && (rmode < (yasymm1+yasymm2)) ) {
4210	imode = 2;
4211	} else if ( (rmode > yasymm1) && (rmode > (yasymm1+yasymm2)) ) {
4212	imode = 3;
4213	}
4214	} while(imode == 0);
4215
4216	// /* determine parameters of the Z distribution */
4217	// force imode (for testing, PK)
4218	// imode = 3;
4219	if (imode == 1) {
4220	z1mean = zheavy1;
4221	z1width = wzasymm1;
4222	}
4223	if (imode == 2) {
4224	z1mean = zheavy2;
4225	z1width = wzasymm2;
4226	}
4227	if (imode == 3) {
4228	z1mean = zsymm;
4229	z1width = wzsymm;
4230	}
4231
4232	if (itest == 1) {
4233	G4cout << "nbre aleatoire tire " << rmode << G4endl;
4234	G4cout << "fission mode " << imode << G4endl;
4235	G4cout << "z1mean= " << z1mean << G4endl;
4236	G4cout << "z1width= " << z1width << G4endl;
4237	}
4238
4239	// /* random decision: Z1 and Z2 at scission: */
4240	z1 = 1.0;
4241	z2 = 1.0;
4242	while ( (z1<5.0) \|\| (z2<5.0) ) {
4243	// Z1 = dble(GAUSSHAZ(K,sngl(Z1mean),sngl(Z1width)));
4244	// z1 = rnd.gaus(z1mean,z1width);
4245	z1 = gausshaz(k, z1mean, z1width);
4246	z2 = z - z1;
4247	}
4248	if (itest == 1) {
4249	G4cout << "ff charge sample " << G4endl;
4250	G4cout << "z1 = " << z1 << G4endl;
4251	G4cout << "z2 = " << z2 << G4endl;
4252	}
4253
4254	// CALL EVEN_ODD(Z1,R_E_O,I_HELP);
4255	// /* Integer proton number with even-odd effect */
4256	// Z1 = REAL(I_HELP)
4257	// /* Z1 = INT(Z1+0.5E0); */
4258	z2 = z - z1;
4259
4260	// /* average N of both fragments: */
4261	if (imode == 1) {
4262	n1mean = (z1 + cpol1 * a/n) * n/z;
4263	}
4264	if (imode == 2) {
4265	n1mean = (z1 + cpol2 * a/n) * n/z;
4266	}
4267	/* CASE(99) ! only for testing;
4268	N1UCD = Z1 * N/Z;
4269	N2UCD = Z2 * N/Z;
4270	re1 = UMASS(Z1,N1UCD,0.6) +;
4271	& UMASS(Z2,N2UCD,0.6) +;
4272	& ECOUL(Z1,N1UCD,0.6,Z2,N2UCD,0.6,d);
4273	re2 = UMASS(Z1,N1UCD+1.,0.6) +;
4274	& UMASS(Z2,N2UCD-1.,0.6) +;
4275	& ECOUL(Z1,N1UCD+1.,0.6,Z2,N2UCD-1.,0.6,d);
4276	re3 = UMASS(Z1,N1UCD+2.,0.6) +;
4277	& UMASS(Z2,N2UCD-2.,0.6) +;
4278	& ECOUL(Z1,N1UCD+2.,0.6,Z2,N2UCD-2.,0.6,d);
4279	eps2 = (re1-2.0*re2+re3) / 2.0;
4280	eps1 = re2 - re1 - eps2;
4281	DN1_POL = - eps1 / (2.0 * eps2);
4282	N1mean = N1UCD + DN1_POL; */
4283	if (imode == 3) {
4284	n1ucd = z1 * n/z;
4285	n2ucd = z2 * n/z;
4286	re1 = umass(z1,n1ucd,0.6) + umass(z2,n2ucd,0.6) + ecoul(z1,n1ucd,0.6,z2,n2ucd,0.6,d);
4287	re2 = umass(z1,n1ucd+1.,0.6) + umass(z2,n2ucd-1.,0.6) + ecoul(z1,n1ucd+1.,0.6,z2,n2ucd-1.,0.6,d);
4288	re3 = umass(z1,n1ucd+2.,0.6) + umass(z2,n2ucd-2.,0.6) + ecoul(z1,n1ucd+2.,0.6,z2,n2ucd-2.,0.6,d);
4289	eps2 = (re1-2.0*re2+re3) / 2.0;
4290	eps1 = re2 - re1 - eps2;
4291	dn1_pol = - eps1 / (2.0 * eps2);
4292	n1mean = n1ucd + dn1_pol;
4293	}
4294	// all fission modes features have been checked CS 11/10/05
4295	n2mean = n - n1mean;
4296	z2mean = z - z1mean;
4297
4298	// /* Excitation energies */
4299	// /* formulated in energies in close consistency with the fission model */
4300
4301	// /* E_defo = UMASS(Z0.5E0,N0.5E0,0.6E0) -
4302	// UMASS(Z0.5E0,N0.5E0,0); */
4303	// /* calculates the deformation energy of the liquid drop for
4304	// deformation beta = 0.6 which is most probable at scission */
4305
4306	// /* N1R and N2R provisionaly taken without fluctuations in
4307	// polarisation: */
4308	n1r = n1mean;
4309	n2r = n2mean;
4310	a1r = n1r + z1;
4311	a2r = n2r + z2;
4312
4313	if (imode == 1) { /* N = 82 */;
4314	//! /* Eexc at scission */
4315	e_scission = max(epsilon_1_scission,1.0);
4316	if (n1mean > (n * 0.5) ) {
4317	//! /* 1. fragment is spherical */
4318	beta1 = 0.0;
4319	beta2 = 0.6;
4320	e1exc = epsilon_1_scission * a1r / a;
4321	e_defo = umass(z2,n2r,beta2) - umass(z2,n2r,0.0);
4322	e2exc = epsilon_1_scission * a2r / a + e_defo;
4323	}
4324	else {
4325	//! /* 2. fragment is spherical */
4326	beta1 = 0.6;
4327	beta2 = 0.0;
4328	e_defo = umass(z1,n1r,beta1) - umass(z1,n1r,0.0);
4329	e1exc = epsilon_1_scission * a1r / a + e_defo;
4330	e2exc = epsilon_1_scission * a2r / a;
4331	}
4332	}
4333
4334	if (imode == 2) {
4335	//! /* N appr. 86 */
4336	e_scission = max(epsilon_2_scission,1.0);
4337	if (n1mean > (n * 0.5) ) {
4338	//! /* 2. fragment is spherical */
4339	beta1 = (n1r - nheavy2) * 0.034 + 0.3;
4340	e_defo = umass(z1,n1r,beta1) - umass(z1,n1r,0.0);
4341	e1exc = epsilon_2_scission * a1r / a + e_defo;
4342	beta2 = 0.6 - beta1;
4343	e_defo = umass(z2,n2r,beta2) - umass(z2,n2r,0.0);
4344	e2exc = epsilon_2_scission * a2r / a + e_defo;
4345	}
4346	else {
4347	//! /* 1. fragment is spherical */
4348	beta2 = (n2r - nheavy2) * 0.034 + 0.3;
4349	e_defo = umass(z2,n2r,beta2) - umass(z2,n2r,0.0);
4350	e2exc = epsilon_2_scission * a2r / a + e_defo;
4351	beta1 = 0.6 - beta2;
4352	e_defo = umass(z1,n1r,beta1) - umass(z1,n1r,0.0);
4353	e1exc = epsilon_2_scission * a1r / a + e_defo;
4354	}
4355	}
4356
4357	if (imode == 3) {
4358	// ! /* Symmetric fission channel */
4359
4360	// /* the fit function for beta is the deformation for
4361	// optimum energy at the scission point, d = 2 */
4362	// /* beta : deformation of symmetric fragments */
4363	// /* beta1 : deformation of first fragment */
4364	// /* beta2 : deformation of second fragment */
4365	beta = 0.177963 + 0.0153241 * zsymm - 0.000162037 * zsymm*zsymm;
4366	beta1 = 0.177963 + 0.0153241 * z1 - 0.000162037 * z1*z1;
4367	// beta1 = 0.6
4368	e_defo1 = umass(z1,n1r,beta1) - umass(z1,n1r,0.0);
4369	beta2 = 0.177963 + 0.0153241 * z2 - 0.000162037 * z2*z2;
4370	// beta2 = 0.6
4371	e_defo2 = umass(z2,n2r,beta2) - umass(z2,n2r,0.0);
4372	e_asym = umass(z1 , n1r, beta1) + umass(z2, n2r ,beta2)
4373	+ ecoul(z1,n1r,beta1,z2,n2r,beta2,2.0)
4374	- 2.0 * umass(zsymm,nsymm,beta)
4375	- ecoul(zsymm,nsymm,beta,zsymm,nsymm,beta,2.0);
4376	// E_asym = CZ_symm * (Z1 - Zsymm)**2
4377	e_scission = max((epsilon_symm_scission - e_asym),1.0);
4378	// /* $LIST(Z1,N1R,Z2,N2R,E_asym,E_scission); */
4379	e1exc = e_scission * a1r / a + e_defo1;
4380	e2exc = e_scission * a2r / a + e_defo2;
4381	}
4382	// Energies checked for all the modes CS 11/10/05
4383
4384	// /* random decision: N1R and N2R at scission, before evaporation: */
4385	// /* CN = UMASS(Zsymm , Nsymm + 1.E0,0) +
4386	// UMASS(Zsymm, Nsymm - 1.E0,0)
4387	// + 1.44E0 * (Zsymm)**2 /
4388	// (r_null*2 ((Asymm+1)1/3 + (Asymm-1)1/3)**2 )
4389	// - 2.E0 * UMASS(Zsymm,Nsymm,0)
4390	// - 1.44E0 * (Zsymm)*2 / (r_null 2.E0 * (Asymm)1/3)2; */
4391
4392
4393	// /* N1width = std::sqrt(0.5E0 * std::sqrt(1.E0/A_levdens(Eld+E_saddle_scission)) / CN); /
4394	// /* 8. 9. 1998: KHS (see also consideration in the first comment block)
4395	// sigma_N(Z=const) = A/Z * sigma_Z(A=const)
4396	// sigma_Z(A=const) = 0.4 to 0.5 (from Lang paper Nucl Phys. A345 (1980) 34)
4397	// sigma_N(Z=const) = 0.45 * A/Z (= 1.16 for 238U)
4398	// therefore: SIGZMIN = 1.16 */
4399
4400	if ( (imode == 1) \|\| (imode == 2) ) {
4401	cn=(umass(z1,n1mean+1.,beta1) + umass(z1,n1mean-1.,beta1)
4402	+ umass(z2,n2mean+1.,beta2) + umass(z2,n2mean-1.,beta2)
4403	+ ecoul(z1,n1mean+1.,beta1,z2,n2mean-1.,beta2,2.0)
4404	+ ecoul(z1,n1mean-1.,beta1,z2,n2mean+1.,beta2,2.0)
4405	- 2.0 * ecoul(z1,n1mean,beta1,z2,n2mean,beta2,2.0)
4406	- 2.0 * umass(z1, n1mean, beta1)
4407	- 2.0 * umass(z2, n2mean, beta2) ) * 0.5;
4408	// /* Coulomb energy neglected for the moment! */
4409	// IF (E_scission.lt.0.) Then
4410	// write(6,*)'<E> E_scission < 0, MODE 1,2'
4411	// ENDIF
4412	// IF (CN.lt.0.) Then
4413	// write(6,*)'CN < 0, MODE 1,2'
4414	// ENDIF
4415	n1width=std::sqrt( (0.5 * (std::sqrt(1.0/a_levdens*(e_scission)))/cn) );
4416	n1width=max(n1width, sigzmin);
4417
4418	// /* random decision: N1R and N2R at scission, before evaporation: */
4419	n1r = 1.0;
4420	n2r = 1.0;
4421	while ( (n1r<5.0) \|\| (n2r<5.0) ) {
4422	// n1r = dble(gausshaz(k,sngl(n1mean),sngl(n1width)));
4423	// n1r = rnd.gaus(n1mean,n1width);
4424	n1r = gausshaz(k, n1mean, n1width);
4425	n2r = n - n1r;
4426	}
4427	// N1R = GAUSSHAZ(K,N1mean,N1width)
4428	if (itest == 1) {
4429	G4cout << "after neutron sample " << n1r << G4endl;
4430	}
4431	n1r = (float)( (int)((n1r+0.5)) );
4432	n2r = n - n1r;
4433
4434	even_odd(z1,r_e_o,i_help);
4435	// /* proton number with even-odd effect */
4436	z1 = (float)(i_help);
4437	z2 = z - z1;
4438
4439	a1r = z1 + n1r;
4440	a2r = z2 + n2r;
4441	}
4442
4443	if (imode == 3) {
4444	//! /* When(3) */
4445	if (nzpol > 0.0) {
4446	// /* We treat a simultaneous split in Z and N to determine polarisation */
4447	cz = ( umass(z1-1., n1mean+1.,beta1)
4448	+ umass(z2+1., n2mean-1.,beta1)
4449	+ umass(z1+1., n1mean-1.,beta2)
4450	+ umass(z2 - 1., n2mean + 1.,beta2)
4451	+ ecoul(z1-1.,n1mean+1.,beta1,z2+1.,n2mean-1.,beta2,2.0)
4452	+ ecoul(z1+1.,n1mean-1.,beta1,z2-1.,n2mean+1.,beta2,2.0)
4453	- 2.0 * ecoul(z1,n1mean,beta1,z2,n2mean,beta2,2.0)
4454	- 2.0 * umass(z1, n1mean,beta1)
4455	- 2.0 * umass(z2, n2mean,beta2) ) * 0.5;
4456	// IF (E_scission.lt.0.) Then
4457	// write(6,*) '<E> E_scission < 0, MODE 1,2'
4458	// ENDIF
4459	// IF (CZ.lt.0.) Then
4460	// write(6,*) 'CZ < 0, MODE 1,2'
4461	// ENDIF
4462	za1width=std::sqrt( (0.5 * std::sqrt(1.0/a_levdens*(e_scission)) / cz) );
4463	za1width=std::sqrt( (max((za1width*za1width-(1.0/12.0)),0.1)) );
4464	// /* Check the value of 0.1 ! */
4465	// /* Shephard correction */
4466	a1r = z1 + n1mean;
4467	a1r = (float)((int)((a1r+0.5)));
4468	a2r = a - a1r;
4469	// /* A1R and A2R are integer numbers now */
4470	// /* $LIST(A1R,A2R,ZA1WIDTH); */
4471
4472	n1ucd = n/a * a1r;
4473	n2ucd = n/a * a2r;
4474	z1ucd = z/a * a1r;
4475	z2ucd = z/a * a2r;
4476
4477	re1 = umass(z1ucd-1.,n1ucd+1.,beta1) + umass(z2ucd+1.,n2ucd-1.,beta2)
4478	+ ecoul(z1ucd-1.,n1ucd+1.,beta1,z2ucd+1.,n2ucd-1.,beta2,d);
4479	re2 = umass(z1ucd,n1ucd,beta1) + umass(z2ucd,n2ucd,beta2)
4480	+ ecoul(z1ucd,n1ucd,beta1,z2ucd,n2ucd,beta2,d);
4481	re3 = umass(z1ucd+1.,n1ucd-1.,beta1) + umass(z2ucd-1.,n2ucd+1.,beta2) +
4482	+ ecoul(z1ucd+1.,n1ucd-1.,beta1,z2ucd-1.,n2ucd+1.,beta2,d);
4483
4484	eps2 = (re1-2.0*re2+re3) / 2.0;
4485	eps1 = (re3 - re1)/2.0;
4486	dn1_pol = - eps1 / (2.0 * eps2);
4487	z1 = z1ucd + dn1_pol;
4488	if (itest == 1) {
4489	G4cout << "before proton sample " << z1 << G4endl;
4490	}
4491	// Z1 = dble(GAUSSHAZ(k,sngl(Z1),sngl(ZA1width)));
4492	// z1 = rnd.gaus(z1,za1width);
4493	z1 = gausshaz(k, z1, za1width);
4494	if (itest == 1) {
4495	G4cout << "after proton sample " << z1 << G4endl;
4496	}
4497	even_odd(z1,r_e_o,i_help);
4498	// /* proton number with even-odd effect */
4499	z1 = (float)(i_help);
4500	z2 = (float)((int)( (z - z1 + 0.5)) );
4501
4502	n1r = a1r - z1;
4503	n2r = n - n1r;
4504	}
4505	else {
4506	// /* First division of protons, then adjustment of neutrons */
4507	cn = ( umass(z1, n1mean+1.,beta1) + umass(z1, n1mean-1., beta1)
4508	+ umass(z2, n2mean+1.,beta2) + umass(z2, n2mean-1., beta2)
4509	+ ecoul(z1,n1mean+1.,beta1,z2,n2mean-1.,beta2,2.0)
4510	+ ecoul(z1,n1mean-1.,beta1,z2,n2mean+1.,beta2,2.0)
4511	- 2.0 * ecoul(z1,n1mean,beta1,z2,n2mean,beta2,2.0)
4512	- 2.0 * umass(z1, n1mean, 0.6)
4513	- 2.0 * umass(z2, n2mean, 0.6) ) * 0.5;
4514	// /* Coulomb energy neglected for the moment! */
4515	// IF (E_scission.lt.0.) Then
4516	// write(6,*) '<E> E_scission < 0, MODE 1,2'
4517	// Endif
4518	// IF (CN.lt.0.) Then
4519	// write(6,*) 'CN < 0, MODE 1,2'
4520	// Endif
4521	n1width=std::sqrt( (0.5 * std::sqrt(1.0/a_levdens*(e_scission)) / cn) );
4522	n1width=max(n1width, sigzmin);
4523
4524	// /* random decision: N1R and N2R at scission, before evaporation: */
4525	// N1R = dble(GAUSSHAZ(k,sngl(N1mean),sngl(N1width)));
4526	// n1r = rnd.gaus(n1mean,n1width);
4527	n1r = gausshaz(k, n1mean, n1width);
4528	n1r = (float)( (int)((n1r+0.5)) );
4529	n2r = n - n1r;
4530
4531	even_odd(z1,r_e_o,i_help);
4532	// /* Integer proton number with even-odd effect */
4533	z1 = (float)(i_help);
4534	z2 = z - z1;
4535
4536	a1r = z1 + n1r;
4537	a2r = z2 + n2r;
4538
4539	}
4540	}
4541
4542	if (itest == 1) {
4543	G4cout << "remid imode = " << imode << G4endl;
4544	G4cout << "n1width = " << n1width << G4endl;
4545	G4cout << "n1r = " << n1r << G4endl;
4546	G4cout << "a1r = " << a1r << G4endl;
4547	G4cout << "n2r = " << n2r << G4endl;
4548	G4cout << "a2r = " << a2r << G4endl;
4549	}
4550	// Up to here: checked CS 11/10/05
4551
4552	// /* Extracted from Lang et al. Nucl. Phys. A 345 (1980) 34 */
4553	e1exc_sigma = 5.5;
4554	e2exc_sigma = 5.5;
4555
4556	neufcentquatrevingtsept:
4557	// E1final = dble(Gausshaz(k,sngl(E1exc),sngl(E1exc_sigma)));
4558	// E2final = dble(Gausshaz(k,sngl(E2exc),sngl(E2exc_sigma)));
4559	// e1final = rnd.gaus(e1exc,e1exc_sigma);
4560	// e2final = rnd.gaus(e2exc,e2exc_sigma);
4561	e1final = gausshaz(k, e1exc, e1exc_sigma);
4562	e2final = gausshaz(k, e2exc, e2exc_sigma);
4563	if ( (e1final < 0.0) \|\| (e2final < 0.0) ) goto neufcentquatrevingtsept;
4564	if (itest == 1) {
4565	G4cout << "sampled exc 1 " << e1final << G4endl;
4566	G4cout << "sampled exc 2 " << e2final << G4endl;
4567	}
4568
4569	// /* OUTPUT QUANTITIES OF THE EVENT GENERATOR: */
4570
4571	// /* Quantities before neutron evaporation */
4572
4573	// /* Neutron number of prefragments: N1R and N2R */
4574	// /* Atomic number of fragments: Z1 and Z2 */
4575	// /* Kinetic energy of fragments: EkinR1, EkinR2 *7
4576
4577	// /* Quantities after neutron evaporation: */
4578
4579	// /* Neutron number of fragments: N1 and N2 */
4580	// /* Mass number of fragments: A1 and A2 */
4581	// /* Atomic number of fragments: Z1 and Z2 */
4582	// /* Number of evaporated neutrons: N1R-N1 and N2R-N2 */
4583	// /* Kinetic energy of fragments: EkinR1*A1/A1R and
4584	// EkinR2A2/A2R /
4585
4586	n1 = n1r;
4587	n2 = n2r;
4588	a1 = n1 + z1;
4589	a2 = n2 + z2;
4590	e1 = e1final;
4591	e2 = e2final;
4592
4593	// /* Pre-neutron-emission total kinetic energy: */
4594	tker = (z1 * z2 * 1.44) /
4595	( r0 * std::pow(a1,0.33333) * (1.0 + 2.0/3.0 * beta1) +
4596	r0 * std::pow(a2,0.33333) * (1.0 + 2.0/3.0 * beta2) + 2.0 );
4597	// /* Pre-neutron-emission kinetic energy of 1. fragment: */
4598	ekinr1 = tker * a2 / a;
4599	// /* Pre-neutron-emission kinetic energy of 2. fragment: */
4600	ekinr2 = tker * a1 / a;
4601
4602	v1 = std::sqrt( (ekinr1/a1) ) * 1.3887;
4603	v2 = std::sqrt( (ekinr2/a2) ) * 1.3887;
4604
4605	if (itest == 1) {
4606	G4cout << "ekinr1 " << ekinr1 << G4endl;
4607	G4cout << "ekinr2 " << ekinr2 << G4endl;
4608	}
4609
4610	milledeux:
4611	//**************************
4612	//*** only symmetric fission
4613	//**************************
4614	// Symmetric fission: Ok! Checked CS 10/10/05
4615	if ( (icz == -1) \|\| (a1 < 0.0) \|\| (a2 < 0.0) ) {
4616	// IF (z.eq.92) THEN
4617	// write(6,*)'symmetric fission'
4618	// write(6,*)'Z,A,E,A1,A2,icz,Atot',Z,A,E,A1,A2,icz,Atot
4619	// END IF
4620
4621	if (itest == 1) {
4622	G4cout << "milledeux: liquid-drop option " << G4endl;
4623	}
4624
4625	n = a-z;
4626	// proton number in symmetric fission (centre) *
4627	zsymm = z / 2.0;
4628	nsymm = n / 2.0;
4629	asymm = nsymm + zsymm;
4630
4631	a_levdens = a / xlevdens;
4632
4633	masscurv = 2.0;
4634	cz_symm = 8.0 / std::pow(z,2) * masscurv;
4635
4636	wzsymm = std::sqrt( (0.5 * std::sqrt(1.0/a_levdens*e) / cz_symm) ) ;
4637
4638	if (itest == 1) {
4639	G4cout << " symmetric high energy fission " << G4endl;
4640	G4cout << "wzsymm " << wzsymm << G4endl;
4641	}
4642
4643	z1mean = zsymm;
4644	z1width = wzsymm;
4645
4646	// random decision: Z1 and Z2 at scission: */
4647	z1 = 1.0;
4648	z2 = 1.0;
4649	while ( (z1 < 5.0) \|\| (z2 < 5.0) ) {
4650	// z1 = dble(gausshaz(kkk,sngl(z1mean),sngl(z1width)));
4651	// z1 = rnd.gaus(z1mean,z1width);
4652	z1 = gausshaz(kkk, z1mean, z1width);
4653	z2 = z - z1;
4654	}
4655
4656	if (itest == 1) {
4657	G4cout << " z1 " << z1 << G4endl;
4658	G4cout << " z2 " << z2 << G4endl;
4659	}
4660	if (itest == 1) {
4661	G4cout << " zsymm " << zsymm << G4endl;
4662	G4cout << " nsymm " << nsymm << G4endl;
4663	G4cout << " asymm " << asymm << G4endl;
4664	}
4665	// CN = UMASS(Zsymm , Nsymm + 1.E0) + UMASS(Zsymm, Nsymm - 1.E0)
4666	// # + 1.44E0 * (Zsymm)**2 /
4667	// # (r_null*2 ((Asymm+1)**(1./3.) +
4668	// # (Asymm-1)(1./3.))2 )
4669	// # - 2.E0 * UMASS(Zsymm,Nsymm)
4670	// # - 1.44E0 * (Zsymm)**2 /
4671	// # (r_null * 2.E0 * (Asymm)(1./3.))2
4672
4673	n1ucd = z1 * n/z;
4674	n2ucd = z2 * n/z;
4675	re1 = umass(z1,n1ucd,0.6) + umass(z2,n2ucd,0.6) +
4676	ecoul(z1,n1ucd,0.6,z2,n2ucd,0.6,2.0);
4677	re2 = umass(z1,n1ucd+1.,0.6) + umass(z2,n2ucd-1.,0.6) +
4678	ecoul(z1,n1ucd+1.,0.6,z2,n2ucd-1.,0.6,2.0);
4679	re3 = umass(z1,n1ucd+2.,0.6) + umass(z2,n2ucd-2.,0.6) +
4680	ecoul(z1,n1ucd+2.,0.6,z2,n2ucd-2.,0.6,2.0);
4681	reps2 = (re1-2.0*re2+re3)/2.0;
4682	reps1 = re2 - re1 -reps2;
4683	rn1_pol = -reps1/(2.0*reps2);
4684	n1mean = n1ucd + rn1_pol;
4685	n2mean = n - n1mean;
4686
4687	if (itest == 1) {
4688	G4cout << " n1mean " << n1mean << G4endl;
4689	G4cout << " n2mean " << n2mean << G4endl;
4690	}
4691
4692	cn = (umass(z1,n1mean+1.,0.0) + umass(z1,n1mean-1.,0.0) +
4693	+ umass(z2,n2mean+1.,0.0) + umass(z2,n2mean-1.,0.0)
4694	- 2.0 * umass(z1,n1mean,0.0) +
4695	- 2.0 * umass(z2,n2mean,0.0) ) * 0.5;
4696	// This is an approximation! Coulomb energy is neglected.
4697
4698	n1width = std::sqrt( (0.5 * std::sqrt(1.0/a_levdens*e) / cn) );
4699
4700	if (itest == 1) {
4701	G4cout << " cn " << cn << G4endl;
4702	G4cout << " n1width " << n1width << G4endl;
4703	}
4704
4705	// random decision: N1R and N2R at scission, before evaporation: */
4706	// N1R = dfloat(NINT(GAUSSHAZ(KKK,sngl(N1mean),sngl(N1width))));
4707	// n1r = (float)( (int)(rnd.gaus(n1mean,n1width)) );
4708	n1r = (float)( (int)(gausshaz(k, n1mean,n1width)) );
4709	n2r = n - n1r;
4710	// Mass of first and second fragment */
4711	a1 = z1 + n1r;
4712	a2 = z2 + n2r;
4713
4714	e1 = e*a1/(a1+a2);
4715	e2 = e - e*a1/(a1+a2);
4716	if (itest == 1) {
4717	G4cout << " n1r " << n1r << G4endl;
4718	G4cout << " n2r " << n2r << G4endl;
4719	}
4720
4721	}
4722
4723	if (itest == 1) {
4724	G4cout << " a1 " << a1 << G4endl;
4725	G4cout << " z1 " << z1 << G4endl;
4726	G4cout << " a2 " << a2 << G4endl;
4727	G4cout << " z2 " << z2 << G4endl;
4728	G4cout << " e1 " << e1 << G4endl;
4729	G4cout << " e2 " << e << G4endl;
4730	}
4731
4732	// /* Pre-neutron-emission total kinetic energy: */
4733	tker = (z1 * z2 * 1.44) /
4734	( r0 * std::pow(a1,0.33333) * (1.0 + 2.0/3.0 * beta1) +
4735	r0 * std::pow(a2,0.33333) * (1.0 + 2.0/3.0 * beta2) + 2.0 );
4736	// /* Pre-neutron-emission kinetic energy of 1. fragment: */
4737	ekin1 = tker * a2 / a;
4738	// /* Pre-neutron-emission kinetic energy of 2. fragment: */
4739	ekin2 = tker * a1 / a;
4740
4741	v1 = std::sqrt( (ekin1/a1) ) * 1.3887;
4742	v2 = std::sqrt( (ekin2/a2) ) * 1.3887;
4743
4744	if (itest == 1) {
4745	G4cout << " kinetic energies " << G4endl;
4746	G4cout << " ekin1 " << ekin1 << G4endl;
4747	G4cout << " ekin2 " << ekin2 << G4endl;
4748	}
4749	}
4750
4751	// SUBROUTINE TRANSLAB(GAMREM,ETREM,CSREM,NOPART,NDEC)
4752	void G4Abla::translab(G4double gamrem, G4double etrem, G4double csrem[4], G4int nopart, G4int ndec)
4753	{
4754	// c Ce subroutine transforme dans un repere 1 les impulsions pcv des
4755	// c particules acv, zcv et de cosinus directeurs xcv, ycv, zcv calculees
4756	// c dans un repere 2.
4757	// c La transformation de lorentz est definie par GAMREM (gamma) et
4758	// c ETREM (eta). La direction du repere 2 dans 1 est donnees par les
4759	// c cosinus directeurs ALREM,BEREM,GAREM (axe oz du repere 2).
4760	// c L'axe oy(2) est fixe par le produit vectoriel oz(1)*oz(2).
4761	// c Le calcul est fait pour les particules de NDEC a iv du common volant.
4762	// C Resultats dans le NTUPLE (common VAR_NTP) decale de NOPART (cascade).
4763
4764	// REAL*8 GAMREM,ETREM,ER,PLABI(3),PLABF(3),R(3,3)
4765	// real*8 MASSE,PTRAV2,CSREM(3),UMA,MELEC,EL
4766	// real*4 acv,zpcv,pcv,xcv,ycv,zcv
4767	// common/volant/acv(200),zpcv(200),pcv(200),xcv(200),
4768	// s ycv(200),zcv(200),iv
4769
4770	// parameter (max=250)
4771	// real*4 EXINI,ENERJ,BIMPACT,PLAB,TETLAB,PHILAB,ESTFIS
4772	// integer AVV,ZVV,JREMN,KFIS,IZFIS,IAFIS
4773	// common/VAR_NTP/MASSINI,MZINI,EXINI,MULNCASC,MULNEVAP,
4774	// +MULNTOT,BIMPACT,JREMN,KFIS,ESTFIS,IZFIS,IAFIS,NTRACK,
4775	// +ITYPCASC(max),AVV(max),ZVV(max),ENERJ(max),PLAB(max),
4776	// +TETLAB(max),PHILAB(max)
4777
4778	// DATA UMA,MELEC/931.4942,0.511/
4779
4780	// C Matrice de rotation dans le labo:
4781	G4double sitet = std::sqrt(std::pow(csrem[1],2)+std::pow(csrem[2],2));
4782	G4double cstet = 0.0, siphi = 0.0, csphi = 0.0;
4783	G4double R[4][4];
4784	for(G4int init_i = 0; init_i < 4; init_i++) {
4785	for(G4int init_j = 0; init_j < 4; init_j++) {
4786	R[init_i][init_j] = 0.0;
4787	}
4788	}
4789
4790	if(sitet > 1.0e-6) { //then
4791	cstet = csrem[3];
4792	siphi = csrem[2]/sitet;
4793	csphi = csrem[1]/sitet;
4794
4795	R[1][1] = cstet*csphi;
4796	R[1][2] = -siphi;
4797	R[1][3] = sitet*csphi;
4798	R[2][1] = cstet*siphi;
4799	R[2][2] = csphi;
4800	R[2][3] = sitet*siphi;
4801	R[3][1] = -sitet;
4802	R[3][2] = 0.0;
4803	R[3][3] = cstet;
4804	}
4805	else {
4806	R[1][1] = 1.0;
4807	R[1][2] = 0.0;
4808	R[1][3] = 0.0;
4809	R[2][1] = 0.0;
4810	R[2][2] = 1.0;
4811	R[2][3] = 0.0;
4812	R[3][1] = 0.0;
4813	R[3][2] = 0.0;
4814	R[3][3] = 1.0;
4815	} //endif
4816
4817	G4int intp = 0;
4818	G4double el = 0.0;
4819	G4double masse = 0.0;
4820	G4double er = 0.0;
4821	G4double plabi[4];
4822	G4double ptrav2 = 0.0;
4823	G4double plabf[4];
4824	G4double bidon = 0.0;
4825	for(G4int init_i = 0; init_i < 4; init_i++) {
4826	plabi[init_i] = 0.0;
4827	plabf[init_i] = 0.0;
4828	}
4829
4830	for(G4int i = ndec; i <= volant->iv; i++) { //do i=ndec,iv
4831	intp = i + nopart;
4832	varntp->ntrack = varntp->ntrack + 1;
4833	if(nint(volant->acv[i]) == 0 && nint(volant->zpcv[i]) == 0) {
4834	if(verboseLevel > 2) {
4835	G4cout <<"Error: Particles with A = 0 Z = 0 detected! " << G4endl;
4836	}
4837	continue;
4838	}
4839	if(varntp->ntrack >= VARNTPSIZE) {
4840	if(verboseLevel > 2) {
4841	G4cout <<"Error! Output data structure not big enough!" << G4endl;
4842	}
4843	}
4844	varntp->avv[intp] = nint(volant->acv[i]);
4845	varntp->zvv[intp] = nint(volant->zpcv[i]);
4846	varntp->itypcasc[intp] = 0;
4847	// transformation de lorentz remnan --> labo:
4848	if (varntp->avv[intp] == -1) { //then
4849	masse = 138.00; // cugnon
4850	// c if (avv(intp).eq.1) masse=938.2796 !cugnon
4851	// c if (avv(intp).eq.4) masse=3727.42 !ok
4852	}
4853	else {
4854	mglms(double(volant->acv[i]),double(volant->zpcv[i]),0, &el);
4855	// assert(isnan(el) == false);
4856	masse = volant->zpcv[i]938.27 + (volant->acv[i] - volant->zpcv[i])939.56 + el;
4857	} //end if
4858
4859	er = std::sqrt(std::pow(volant->pcv[i],2) + std::pow(masse,2));
4860	// assert(isnan(er) == false);
4861	plabi[1] = volant->pcv[i]*(volant->xcv[i]);
4862	plabi[2] = volant->pcv[i]*(volant->ycv[i]);
4863	plabi[3] = eretrem + gamrem(volant->pcv[i])*(volant->zcv[i]);
4864
4865	ptrav2 = std::pow(plabi[1],2) + std::pow(plabi[2],2) + std::pow(plabi[3],2);
4866	// assert(isnan(ptrav2) == false);
4867	varntp->plab[intp] = std::sqrt(ptrav2);
4868	varntp->enerj[intp] = std::sqrt(ptrav2 + std::pow(masse,2)) - masse;
4869
4870	// Rotation dans le labo:
4871	for(G4int j = 1; j <= 3; j++) { //do j=1,3
4872	plabf[j] = 0.0;
4873	for(G4int k = 1; k <= 3; k++) { //do k=1,3
4874	plabf[j] = plabf[j] + R[k][j]*plabi[k]; // :::Fixme::: (indices?)
4875	} // end do
4876	} // end do
4877	// C impulsions dans le nouveau systeme copiees dans /volant/
4878	volant->pcv[i] = varntp->plab[intp];
4879	ptrav2 = std::sqrt(std::pow(plabf[1],2) + std::pow(plabf[2],2) + std::pow(plabf[3],2));
4880	if(ptrav2 >= 1.0e-6) { //then
4881	volant->xcv[i] = plabf[1]/ptrav2;
4882	volant->ycv[i] = plabf[2]/ptrav2;
4883	volant->zcv[i] = plabf[3]/ptrav2;
4884	}
4885	else {
4886	volant->xcv[i] = 1.0;
4887	volant->ycv[i] = 0.0;
4888	volant->zcv[i] = 0.0;
4889	} //endif
4890	// impulsions dans le nouveau systeme copiees dans /VAR_NTP/
4891	if(varntp->plab[intp] >= 1.0e-6) { //then
4892	bidon = plabf[3]/(varntp->plab[intp]);
4893	// assert(isnan(bidon) == false);
4894	if(bidon > 1.0) {
4895	bidon = 1.0;
4896	}
4897	if(bidon < -1.0) {
4898	bidon = -1.0;
4899	}
4900	varntp->tetlab[intp] = std::acos(bidon);
4901	sitet = std::sin(varntp->tetlab[intp]);
4902	varntp->philab[intp] = std::atan2(plabf[2],plabf[1]);
4903	varntp->tetlab[intp] = varntp->tetlab[intp]*57.2957795;
4904	varntp->philab[intp] = varntp->philab[intp]*57.2957795;
4905	}
4906	else {
4907	varntp->tetlab[intp] = 90.0;
4908	varntp->philab[intp] = 0.0;
4909	} // endif
4910	} // end do
4911	}
4912	// C-------------------------------------------------------------------------
4913
4914	// SUBROUTINE TRANSLABPF(MASSE1,T1,P1,CTET1,PHI1,GAMREM,ETREM,R,
4915	// s PLAB1,GAM1,ETA1,CSDIR)
4916	void G4Abla::translabpf(G4double masse1, G4double t1, G4double p1, G4double ctet1,
4917	G4double phi1, G4double gamrem, G4double etrem, G4double R[][4],
4918	G4double plab1, G4double gam1, G4double *eta1, G4double csdir[])
4919	{
4920	// C Calcul de l'impulsion du PF (PLAB1, cos directeurs CSDIR(3)) dans le
4921	// C systeme remnant et des coefs de Lorentz GAM1,ETA1 de passage
4922	// c du systeme PF --> systeme remnant.
4923	// c
4924	// C Input: MASSE1, T1 (energie cinetique), CTET1,PHI1 (cosTHETA et PHI)
4925	// C (le PF dans le systeme du Noyau de Fission (NF)).
4926	// C GAMREM,ETREM les coefs de Lorentz systeme NF --> syst remnant,
4927	// C R(3,3) la matrice de rotation systeme NF--> systeme remnant.
4928	// C
4929	// C
4930	// REAL*8 MASSE1,T1,P1,CTET1,PHI1,GAMREM,ETREM,R(3,3),
4931	// s PLAB1,GAM1,ETA1,CSDIR(3),ER,SITET,PLABI(3),PLABF(3)
4932
4933	G4double er = t1 + masse1;
4934
4935	G4double sitet = std::sqrt(1.0 - std::pow(ctet1,2));
4936
4937	G4double plabi[4];
4938	G4double plabf[4];
4939	for(G4int init_i = 0; init_i < 4; init_i++) {
4940	plabi[init_i] = 0.0;
4941	plabf[init_i] = 0.0;
4942	}
4943
4944	// C ----Transformation de Lorentz Noyau fissionnant --> Remnant:
4945	plabi[1] = p1sitetstd::cos(phi1);
4946	plabi[2] = p1sitetstd::sin(phi1);
4947	plabi[3] = eretrem + gamremp1*ctet1;
4948
4949	// C ----Rotation du syst Noyaut Fissionant vers syst remnant:
4950	for(G4int j = 1; j <= 3; j++) { // do j=1,3
4951	plabf[j] = 0.0;
4952	for(G4int k = 1; k <= 3; k++) { //do k=1,3
4953	// plabf[j] = plabf[j] + R[j][k]*plabi[k];
4954	plabf[j] = plabf[j] + R[k][j]*plabi[k];
4955	} //end do
4956	} //end do
4957	// C ----Cosinus directeurs et coefs de la transf de Lorentz dans le
4958	// c nouveau systeme:
4959	(*plab1) = std::pow(plabf[1],2) + std::pow(plabf[2],2) + std::pow(plabf[3],2);
4960	(gam1) = std::sqrt(std::pow(masse1,2) + (plab1))/masse1;
4961	(plab1) = std::sqrt((plab1));
4962	(eta1) = (plab1)/masse1;
4963
4964	if((*plab1) <= 1.0e-6) { //then
4965	csdir[1] = 0.0;
4966	csdir[2] = 0.0;
4967	csdir[3] = 1.0;
4968	}
4969	else {
4970	for(G4int i = 1; i <= 3; i++) { //do i=1,3
4971	csdir[i] = plabf[i]/(*plab1);
4972	} // end do
4973	} //endif
4974	}
4975
4976	// SUBROUTINE LOR_AB(GAM,ETA,Ein,Pin,Eout,Pout)
4977	void G4Abla::lorab(G4double gam, G4double eta, G4double ein, G4double pin[],
4978	G4double *eout, G4double pout[])
4979	{
4980	// C Transformation de lorentz brute pour vÃ©rifs.
4981	// C P(3) = P_longitudinal (transformÃ©)
4982	// C P(1) et P(2) = P_transvers (non transformÃ©s)
4983	// DIMENSION Pin(3),Pout(3)
4984	// REAL*8 GAM,ETA,Ein
4985
4986	pout[1] = pin[1];
4987	pout[2] = pin[2];
4988	(eout) = gamein + eta*pin[3];
4989	pout[3] = etaein + gampin[3];
4990	}
4991
4992	// SUBROUTINE ROT_AB(R,Pin,Pout)
4993	void G4Abla::rotab(G4double R[4][4], G4double pin[4], G4double pout[4])
4994	{
4995	// C Rotation d'un vecteur
4996	// DIMENSION Pin(3),Pout(3)
4997	// REAL*8 R(3,3)
4998
4999	for(G4int i = 1; i <= 3; i++) { // do i=1,3
5000	pout[i] = 0.0;
5001	for(G4int j = 1; j <= 3; j++) { //do j=1,3
5002	// pout[i] = pout[i] + R[i][j]*pin[j];
5003	pout[i] = pout[i] + R[j][i]*pin[j];
5004	} // enddo
5005	} //enddo
5006	}
5007
5008	// Methods related to the internal ABLA random number generator. In
5009	// the future the random number generation must be factored into its
5010	// own class
5011
5012	void G4Abla::standardRandom(G4double rndm, G4long seed)
5013	{
5014	(seed) = (seed); // Avoid warning during compilation.
5015	// Use Geant4 G4UniformRand
5016	(*rndm) = G4UniformRand();
5017	}
5018
5019	G4double G4Abla::haz(G4int k)
5020	{
5021	const G4int pSize = 110;
5022	static G4double p[pSize];
5023	static G4long ix = 0, i = 0;
5024	static G4double x = 0.0, y = 0.0, a = 0.0, haz = 0.0;
5025	// k =< -1 on initialise
5026	// k = -1 c'est reproductible
5027	// k < -1 \|\| k > -1 ce n'est pas reproductible
5028
5029	// Zero is invalid random seed. Set proper value from our random seed collection:
5030	if(ix == 0) {
5031	ix = hazard->ial;
5032	}
5033
5034	if (k <= -1) { //then
5035	if(k == -1) { //then
5036	ix = 0;
5037	}
5038	else {
5039	x = 0.0;
5040	y = secnds(int(x));
5041	ix = int(y * 100 + 43543000);
5042	if(mod(ix,2) == 0) {
5043	ix = ix + 1;
5044	}
5045	}
5046
5047	// Here we are using random number generator copied from INCL code
5048	// instead of the CERNLIB one! This causes difficulties for
5049	// automatic testing since the random number generators, and thus
5050	// the behavior of the routines in C++ and FORTRAN versions is no
5051	// longer exactly the same!
5052	standardRandom(&x, &ix);
5053	for(G4int i = 0; i < pSize; i++) { //do i=1,110
5054	standardRandom(&(p[i]), &ix);
5055	}
5056	standardRandom(&a, &ix);
5057	k = 0;
5058	}
5059
5060	i = nint(100*a)+1;
5061	haz = p[i];
5062	standardRandom(&a, &ix);
5063	p[i] = a;
5064
5065	hazard->ial = ix;
5066
5067	return haz;
5068	}
5069
5070
5071	G4double G4Abla::gausshaz(int k, double xmoy, double sig)
5072	{
5073	// Gaussian random numbers:
5074
5075	// 1005 C*** TIRAGE ALEATOIRE DANS UNE GAUSSIENNE DE LARGEUR SIG ET MOYENNE XMOY
5076	static G4int iset = 0;
5077	static G4double v1,v2,r,fac,gset,gausshaz;
5078
5079	if(iset == 0) { //then
5080	do {
5081	v1 = 2.0*haz(k) - 1.0;
5082	v2 = 2.0*haz(k) - 1.0;
5083	r = std::pow(v1,2) + std::pow(v2,2);
5084	} while(r >= 1);
5085
5086	fac = std::sqrt(-2.*std::log(r)/r);
5087	// assert(isnan(fac) == false);
5088	gset = v1*fac;
5089	gausshaz = v2facsig+xmoy;
5090	iset = 1;
5091	}
5092	else {
5093	gausshaz=gset*sig+xmoy;
5094	iset=0;
5095	}
5096
5097	return gausshaz;
5098	}
5099
5100
5101	// Utilities
5102
5103	G4double G4Abla::min(G4double a, G4double b)
5104	{
5105	if(a < b) {
5106	return a;
5107	}
5108	else {
5109	return b;
5110	}
5111	}
5112
5113	G4int G4Abla::min(G4int a, G4int b)
5114	{
5115	if(a < b) {
5116	return a;
5117	}
5118	else {
5119	return b;
5120	}
5121	}
5122
5123	G4double G4Abla::max(G4double a, G4double b)
5124	{
5125	if(a > b) {
5126	return a;
5127	}
5128	else {
5129	return b;
5130	}
5131	}
5132
5133	G4int G4Abla::max(G4int a, G4int b)
5134	{
5135	if(a > b) {
5136	return a;
5137	}
5138	else {
5139	return b;
5140	}
5141	}
5142
5143	G4int G4Abla::nint(G4double number)
5144	{
5145	G4double intpart = 0.0;
5146	G4double fractpart = 0.0;
5147	fractpart = std::modf(number, &intpart);
5148	if(number == 0) {
5149	return 0;
5150	}
5151	if(number > 0) {
5152	if(fractpart < 0.5) {
5153	return int(std::floor(number));
5154	}
5155	else {
5156	return int(std::ceil(number));
5157	}
5158	}
5159	if(number < 0) {
5160	if(fractpart < -0.5) {
5161	return int(std::floor(number));
5162	}
5163	else {
5164	return int(std::ceil(number));
5165	}
5166	}
5167
5168	return int(std::floor(number));
5169	}
5170
5171	G4int G4Abla::secnds(G4int x)
5172	{
5173	time_t mytime;
5174	tm *mylocaltime;
5175
5176	time(&mytime);
5177	mylocaltime = localtime(&mytime);
5178
5179	if(x == 0) {
5180	return(mylocaltime->tm_hour6060 + mylocaltime->tm_min*60 + mylocaltime->tm_sec);
5181	}
5182	else {
5183	return(mytime - x);
5184	}
5185	}
5186
5187	G4int G4Abla::mod(G4int a, G4int b)
5188	{
5189	if(b != 0) {
5190	return (a - (a/b)*b);
5191	}
5192	else {
5193	return 0;
5194	}
5195	}
5196
5197	G4double G4Abla::dmod(G4double a, G4double b)
5198	{
5199	if(b != 0) {
5200	return (a - (a/b)*b);
5201	}
5202	else {
5203	return 0.0;
5204	}
5205	}
5206
5207	G4double G4Abla::dint(G4double a)
5208	{
5209	G4double value = 0.0;
5210
5211	if(a < 0.0) {
5212	value = double(std::ceil(a));
5213	}
5214	else {
5215	value = double(std::floor(a));
5216	}
5217
5218	return value;
5219	}
5220
5221	G4int G4Abla::idint(G4double a)
5222	{
5223	G4int value = 0;
5224
5225	if(a < 0) {
5226	value = int(std::ceil(a));
5227	}
5228	else {
5229	value = int(std::floor(a));
5230	}
5231
5232	return value;
5233	}
5234
5235	G4int G4Abla::idnint(G4double value)
5236	{
5237	G4double valueCeil = int(std::ceil(value));
5238	G4double valueFloor = int(std::floor(value));
5239
5240	if(std::fabs(value - valueCeil) < std::fabs(value - valueFloor)) {
5241	return int(valueCeil);
5242	}
5243	else {
5244	return int(valueFloor);
5245	}
5246	}
5247
5248	G4double G4Abla::dmin1(G4double a, G4double b, G4double c)
5249	{
5250	if(a < b && a < c) {
5251	return a;
5252	}
5253	if(b < a && b < c) {
5254	return b;
5255	}
5256	if(c < a && c < b) {
5257	return c;
5258	}
5259	return a;
5260	}
5261
5262	G4double G4Abla::utilabs(G4double a)
5263	{
5264	if(a > 0) {
5265	return a;
5266	}
5267	if(a < 0) {
5268	return (-1*a);
5269	}
5270	if(a == 0) {
5271	return a;
5272	}
5273
5274	return a;
5275	}

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

Context Navigation

source: trunk/source/processes/hadronic/models/abla/src/G4Abla.cc @ 1350

Download in other formats: