Context Navigation

source: trunk/source/processes/hadronic/models/incl/src/G4Abla.cc @ 847

Last change on this file since 847 was 819, checked in by garnier, 16 years ago
import all except CVS
File size: 188.7 KB

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

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

Download in other formats: