Context Navigation

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

Last change on this file since 1239 was 1196, checked in by garnier, 15 years ago
update CVS release candidate geant4.9.3.01
File size: 140.2 KB

Line
1	//
2	// ********************************************************************
3	// * License and Disclaimer *
4	// * *
5	// * The Geant4 software is copyright of the Copyright Holders of *
6	// * the Geant4 Collaboration. It is provided under the terms and *
7	// * conditions of the Geant4 Software License, included in the file *
8	// * LICENSE and available at http://cern.ch/geant4/license . These *
9	// * include a list of copyright holders. *
10	// * *
11	// * Neither the authors of this software system, nor their employing *
12	// * institutes,nor the agencies providing financial support for this *
13	// * work make any representation or warranty, express or implied, *
14	// * regarding this software system or assume any liability for its *
15	// * use. Please see the license in the file LICENSE and URL above *
16	// * for the full disclaimer and the limitation of liability. *
17	// * *
18	// * This code implementation is the result of the scientific and *
19	// * technical work of the GEANT4 collaboration. *
20	// * By using, copying, modifying or distributing the software (or *
21	// * any work based on the software) you agree to acknowledge its *
22	// * use in resulting scientific publications, and indicate your *
23	// * acceptance of all terms of the Geant4 Software license. *
24	// ********************************************************************
25	//
26	// $Id: G4Abla.cc,v 1.20 2009/11/18 10:43:14 kaitanie 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 "G4InclRandomNumbers.hh"
39	#include "G4Ranecu.hh"
40	#include "G4AblaFissionSimfis18.hh"
41	#include "G4AblaFission.hh"
42
43	G4Abla::G4Abla()
44	{
45	verboseLevel = 0;
46	ilast = 0;
47	}
48
49	G4Abla::G4Abla(G4Hazard hazard, G4Volant volant)
50	{
51	verboseLevel = 0;
52	ilast = 0;
53	volant = volant; // ABLA internal particle data
54	volant->iv = 0;
55	hazard = hazard; // Random seeds
56
57	randomGenerator = new G4InclGeant4Random();
58	//randomGenerator = new G4Ranecu();
59	varntp = new G4VarNtp();
60	pace = new G4Pace();
61	ald = new G4Ald();
62	eenuc = new G4Eenuc();
63	ec2sub = new G4Ec2sub();
64	ecld = new G4Ecld();
65	fb = new G4Fb();
66	fiss = new G4Fiss();
67	opt = new G4Opt();
68	}
69
70	G4Abla::G4Abla(G4Hazard aHazard, G4Volant aVolant, G4VarNtp *aVarntp)
71	{
72	verboseLevel = 0;
73	ilast = 0;
74	volant = aVolant; // ABLA internal particle data
75	volant->iv = 0;
76	hazard = aHazard; // Random seeds
77	varntp = aVarntp; // Output data structure
78	varntp->ntrack = 0;
79
80	randomGenerator = new G4InclGeant4Random();
81	//randomGenerator = new G4Ranecu();
82
83	// ABLA fission
84	// fissionModel = new G4AblaFissionSimfis18(hazard, randomGenerator);
85	fissionModel = new G4AblaFission(hazard, randomGenerator);
86	if(verboseLevel > 0) {
87	fissionModel->about();
88	}
89	verboseLevel = 0;
90	pace = new G4Pace();
91	ald = new G4Ald();
92	eenuc = new G4Eenuc();
93	ec2sub = new G4Ec2sub();
94	ecld = new G4Ecld();
95	fb = new G4Fb();
96	fiss = new G4Fiss();
97	opt = new G4Opt();
98	}
99
100	G4Abla::~G4Abla()
101	{
102	delete randomGenerator;
103	delete pace;
104	delete ald;
105	delete eenuc;
106	delete ec2sub;
107	delete ecld;
108	delete fb;
109	delete fiss;
110	delete opt;
111	}
112
113	// Main interface to the evaporation
114
115	// Possible problem with generic Geant4 interface: ABLA evaporation
116	// needs angular momentum information (calculated by INCL) to
117	// work. Maybe there is a way to obtain this information from
118	// G4Fragment?
119
120	void G4Abla::breakItUp(G4double nucleusA, G4double nucleusZ, G4double nucleusMass, G4double excitationEnergy,
121	G4double angularMomentum, G4double recoilEnergy, G4double momX, G4double momY, G4double momZ,
122	G4int eventnumber)
123	{
124	const G4double uma = 931.4942;
125	const G4double melec = 0.511;
126	const G4double fmp = 938.27231;
127	const G4double fmn = 939.56563;
128
129	G4double alrem = 0.0, berem = 0.0, garem = 0.0;
130	G4double R[4][4]; // Rotation matrix
131	G4double csdir1[4];
132	G4double csdir2[4];
133	G4double csrem[4];
134	G4double pfis_rem[4];
135	G4double pf1_rem[4];
136	for(G4int init_i = 0; init_i < 4; init_i++) {
137	csdir1[init_i] = 0.0;
138	csdir2[init_i] = 0.0;
139	csrem[init_i] = 0.0;
140	pfis_rem[init_i] = 0.0;
141	pf1_rem[init_i] = 0.0;
142	for(G4int init_j = 0; init_j < 4; init_j++) {
143	R[init_i][init_j] = 0.0;
144	}
145	}
146
147	G4double plab1 = 0.0, gam1 = 0.0, eta1 = 0.0;
148	G4double plab2 = 0.0, gam2 = 0.0, eta2 = 0.0;
149
150	G4double sitet = 0.0;
151	G4double stet1 = 0.0;
152	G4double stet2 = 0.0;
153
154	G4int nbpevap = 0;
155	G4int mempaw = 0, memiv = 0;
156
157	G4double e_evapo = 0.0;
158	G4double el = 0.0;
159	G4double fmcv = 0.0;
160
161	G4double aff1 = 0.0;
162	G4double zff1 = 0.0;
163	G4double eff1 = 0.0;
164	G4double aff2 = 0.0;
165	G4double zff2 = 0.0;
166	G4double eff2 = 0.0;
167
168	G4double v1 = 0.0, v2 = 0.0;
169
170	G4double t2 = 0.0;
171	G4double ctet1 = 0.0;
172	G4double ctet2 = 0.0;
173	G4double phi1 = 0.0;
174	G4double phi2 = 0.0;
175	G4double p2 = 0.0;
176	G4double epf2_out = 0.0 ;
177	G4int lma_pf1 = 0, lmi_pf1 = 0;
178	G4int lma_pf2 = 0, lmi_pf2 = 0;
179	G4int nopart = 0;
180
181	G4double cst = 0.0, sst = 0.0, csf = 0.0, ssf = 0.0;
182
183	G4double zf = 0.0, af = 0.0, mtota = 0.0, pleva = 0.0, pxeva = 0.0, pyeva = 0.0;
184	G4int ff = 0;
185	G4int inum = eventnumber;
186	G4int inttype = 0;
187	G4double esrem = excitationEnergy;
188
189	G4double aprf = nucleusA;
190	G4double zprf = nucleusZ;
191	G4double mcorem = nucleusMass;
192	G4double ee = excitationEnergy;
193	G4double jprf = angularMomentum; // actually root-mean-squared
194
195	G4double erecrem = recoilEnergy;
196	G4double trem = 0.0;
197	G4double pxrem = momX;
198	G4double pyrem = momY;
199	G4double pzrem = momZ;
200
201	G4double remmass = 0.0;
202
203	volant->clear(); // Clean up an initialize ABLA output.
204	varntp->clear(); // Clean up an initialize ABLA output.
205	varntp->ntrack = 0;
206	volant->iv = 0;
207	//volant->iv = 1;
208
209	G4double pcorem = std::sqrt(erecrem(erecrem +2.938.2796*nucleusA));
210	// G4double pcorem = std::sqrt(std::pow(momX,2) + std::pow(momY,2) + std::pow(momZ,2));
211	if(pcorem != 0) { // Guard against division by zero.
212	alrem = pxrem/pcorem;
213	berem = pyrem/pcorem;
214	garem = pzrem/pcorem;
215	} else {
216	alrem = 0.0;
217	berem = 0.0;
218	garem = 0.0;
219	}
220
221	G4int idebug = 0;
222	G4int itest = 0;
223	if(idebug == 1) {
224	zprf = 81.;
225	aprf = 201.;
226	// ee = 86.5877686;
227	ee = 300.0;
228	jprf = 32.;
229	zf = 0.;
230	af = 0.;
231	mtota = 0.;
232	pleva = 0.;
233	pxeva = 0.;
234	pyeva = 0.;
235	ff = -1;
236	inttype = 0;
237	inum = 2;
238	}
239	if(itest == 1) {
240	G4cout <<" PK::: EVAPORA event " << inum << G4endl;
241	G4cout <<" PK::: zprf = " << zprf << G4endl;
242	G4cout <<" PK::: aprf = " << aprf << G4endl;
243	G4cout <<" PK::: ee = " << ee << G4endl;
244	G4cout <<" PK::: eeDiff = " << ee - 86.5877686 << G4endl;
245	G4cout <<" PK::: jprf = " << jprf << G4endl;
246	G4cout <<" PK::: zf = " << zf << G4endl;
247	G4cout <<" PK::: af = " << af << G4endl;
248	G4cout <<" PK::: mtota = " << mtota << G4endl;
249	G4cout <<" PK::: pleva = " << pleva << G4endl;
250	G4cout <<" PK::: pxeva = " << pxeva << G4endl;
251	G4cout <<" PK::: pyeva = " << pyeva << G4endl;
252	G4cout <<" PK::: ff = " << ff << G4endl;
253	G4cout <<" PK::: inttype = " << inttype << G4endl;
254	G4cout <<" PK::: inum = " << inum << G4endl;
255	}
256
257	if(esrem >= 1.0e-3) {
258	evapora(zprf,aprf,ee,jprf, &zf, &af, &mtota, &pleva, &pxeva, &pyeva, &ff, &inttype, &inum);
259	}
260	else {
261	ff = 0;
262	zf = zprf;
263	af = aprf;
264	pxeva = pxrem;
265	pyeva = pyrem;
266	pleva = pzrem;
267	}
268
269	if (ff == 1) {
270	// Fission:
271	// variable ee: Energy of fissioning nucleus above the fission barrier.
272	// Calcul des impulsions des particules evaporees (avant fission)
273	// dans le systeme labo.
274
275	if(verboseLevel > 2)
276	G4cout << __FILE__ << ":" << __LINE__ << " Entering fission code." << G4endl;
277	trem = double(erecrem);
278	remmass = pace2(aprf,zprf) + aprfuma - zprfmelec; // canonic
279	// remmass = mcorem + double(esrem); // ok
280	// remmass = mcorem; //cugnon
281	varntp->kfis = 1;
282	G4double gamrem = (remmass + trem)/remmass;
283	G4double etrem = std::sqrt(trem(trem + 2.0remmass))/remmass;
284
285	// This is not treated as accurately as for the non fission case for which
286	// the remnant mass is computed to satisfy the energy conservation
287	// of evaporated particles. But it is not bad and more canonical!
288	remmass = pace2(aprf,zprf) + aprfuma - zprfmelec+double(esrem); // !canonic
289	// Essais avec la masse de KHS (9/2002):
290	el = 0.0;
291	mglms(aprf,zprf,0,&el);
292	remmass = zprffmp + (aprf-zprf)fmn + el + double(esrem);
293	gamrem = std::sqrt(std::pow(pcorem,2) + std::pow(remmass,2))/remmass;
294	etrem = pcorem/remmass;
295
296	csrem[0] = 0.0; // Should not be used.
297	csrem[1] = alrem;
298	csrem[2] = berem;
299	csrem[3] = garem;
300	if(verboseLevel > 1) {
301	G4cout << __FILE__ << ":" << __LINE__ << " momRem = (" << pxrem << ", " << pyrem << ", " << pzrem << ")" << G4endl;
302	G4cout << __FILE__ << ":" << __LINE__ << " csrem = (" << csrem[1] << ", " << csrem[2] << ", " << csrem[3] << ")" << G4endl;
303	}
304
305	// C Pour VÃ©rif Remnant = evapo(Pre fission) + Noyau_fissionant (systÃšme Remnant)
306	G4double bil_e = 0.0;
307	G4double bil_px = 0.0;
308	G4double bil_py = 0.0;
309	G4double bil_pz = 0.0;
310	G4double masse = 0.0;
311
312	for(G4int iloc = 1; iloc <= volant->iv; iloc++) { //DO iloc=1,iv
313	mglms(double(volant->acv[iloc]),double(volant->zpcv[iloc]),0,&el);
314	masse = volant->zpcv[iloc]fmp + (volant->acv[iloc] - volant->zpcv[iloc])fmn + el;
315	bil_e = bil_e + std::sqrt(std::pow(volant->pcv[iloc],2) + std::pow(masse,2));
316	bil_px = bil_px + volant->pcv[iloc]*(volant->xcv[iloc]);
317	bil_py = bil_py + volant->pcv[iloc]*(volant->ycv[iloc]);
318	bil_pz = bil_pz + volant->pcv[iloc]*(volant->zcv[iloc]);
319	} // enddo
320	// C Ce bilan (impulsion nulle) est parfait. (Bil_Px=Bil_Px+PXEVA....)
321
322	G4int ndec = 1;
323
324	if(volant->iv != 0) { //then
325	if(verboseLevel > 2) {
326	G4cout <<"varntp->ntrack = " << varntp->ntrack << G4endl;
327	G4cout <<"1st Translab: Adding indices from " << ndec << " to " << volant->iv << G4endl;
328	}
329	nopart = varntp->ntrack - 1;
330	translab(gamrem,etrem,csrem,nopart,ndec);
331	if(verboseLevel > 2) {
332	G4cout <<"Translab complete!" << G4endl;
333	G4cout <<"varntp->ntrack = " << varntp->ntrack << G4endl;
334	}
335	}
336	nbpevap = volant->iv; // nombre de particules d'evaporation traitees
337
338	// C
339	// C Now calculation of the fission fragment distribution including
340	// C evaporation from the fragments.
341	// C
342
343	// C Distribution of the fission fragments:
344
345	// void fissionDistri(G4double a,G4double z,G4double e,
346	// G4double &a1,G4double &z1,G4double &e1,G4double &v1,
347	// G4double &a2,G4double &z2,G4double &e2,G4double &v2);
348
349	//fissionModel->fissionDistri(af,zf,ee,aff1,zff1,eff1,v1,aff2,zff2,eff2,v2);
350	fissionModel->doFission(af,zf,ee,aff1,zff1,eff1,v1,aff2,zff2,eff2,v2);
351	// C verif des A et Z decimaux:
352	G4int na_f = int(std::floor(af + 0.5));
353	G4int nz_f = int(std::floor(zf + 0.5));
354	varntp->izfis = nz_f; // copie dans le ntuple
355	varntp->iafis = na_f;
356	G4int na_pf1 = int(std::floor(aff1 + 0.5));
357	G4int nz_pf1 = int(std::floor(zff1 + 0.5));
358	G4int na_pf2 = int(std::floor(aff2 + 0.5));
359	G4int nz_pf2 = int(std::floor(zff2 + 0.5));
360
361	if((na_f != (na_pf1+na_pf2)) \|\| (nz_f != (nz_pf1+nz_pf2))) {
362	if(verboseLevel > 2) {
363	G4cout <<"problemes arrondis dans la fission " << G4endl;
364	G4cout << "af,zf,aff1,zff1,aff2,zff2" << G4endl;
365	G4cout << af <<" , " << zf <<" , " << aff1 <<" , " << zff1 <<" , " << aff2 <<" , " << zff2 << G4endl;
366	G4cout << "a,z,a1,z1,a2,z2 integer" << G4endl;
367	G4cout << na_f <<" , " << nz_f <<" , " << na_pf1 <<" , " << nz_pf1 <<" , " << na_pf2 <<" , " << nz_pf2 << G4endl;
368	}
369	}
370
371	// Calcul de l'impulsion des PF dans le syteme noyau de fission:
372	G4int kboud = idnint(zf);
373	G4int jboud = idnint(af-zf);
374	//G4double ef = fb->efa[kboud][jboud]; // barriere de fission
375	G4double ef = fb->efa[jboud][kboud]; // barriere de fission
376	varntp->estfis = ee + ef; // copie dans le ntuple
377
378	// C MASSEF = pace2(AF,ZF)
379	// C MASSEF = MASSEF + AFUMA - ZFMELEC + EE + EF
380	// C MASSE1 = pace2(DBLE(AFF1),DBLE(ZFF1))
381	// C MASSE1 = MASSE1 + AFF1UMA - ZFF1MELEC + EFF1
382	// C MASSE2 = pace2(DBLE(AFF2),DBLE(ZFF2))
383	// C MASSE2 = MASSE2 + AFF2UMA - ZFF2MELEC + EFF2
384	// C WRITE(6,*) 'MASSEF,MASSE1,MASSE2',MASSEF,MASSE1,MASSE2
385	// C MGLMS est la fonction de masse cohÃ©rente avec KHS evapo-fis.
386	// C Attention aux parametres, ici 0=OPTSHP, NO microscopic correct.
387	mglms(af,zf,0,&el);
388	G4double massef = zffmp + (af - zf)fmn + el + ee + ef;
389	mglms(double(aff1),double(zff1),0,&el);
390	G4double masse1 = zff1fmp + (aff1-zff1)fmn + el + eff1;
391	mglms(aff2,zff2,0,&el);
392	G4double masse2 = zff2fmp + (aff2 - zff2)fmn + el + eff2;
393	// C WRITE(6,*) 'MASSEF,MASSE1,MASSE2',MASSEF,MASSE1,MASSE2
394	G4double b = massef - masse1 - masse2;
395	if(b < 0.0) { //then
396	b=0.0;
397	if(verboseLevel > 2) {
398	G4cout <<"anomalie dans la fission: " << G4endl;
399	G4cout << inum<< " , " << af<< " , " <<zf<< " , " <<massef<< " , " <<aff1<< " , " <<zff1<< " , " <<masse1<< " , " <<aff2<< " , " <<zff2<< " , " << masse2 << G4endl;
400	}
401	} //endif
402	G4double t1 = b(b + 2.0masse2)/(2.0*massef);
403	G4double p1 = std::sqrt(t1(t1 + 2.0masse1));
404
405	G4double rndm;
406	rndm = randomGenerator->getRandom();
407	// standardRandom(&rndm, &(hazard->igraine[13]));
408	ctet1 = 2.0*rndm - 1.0;
409	// standardRandom(&rndm,&(hazard->igraine[9]));
410	rndm = randomGenerator->getRandom();
411	phi1 = rndm2.03.141592654;
412
413	// C ----Coefs de la transformation de Lorentz (noyau de fission -> Remnant)
414	G4double peva = std::pow(pxeva,2) + std::pow(pyeva,2) + std::pow(pleva,2);
415	G4double gamfis = std::sqrt(std::pow(massef,2) + peva)/massef;
416	peva = std::sqrt(peva);
417	G4double etfis = peva/massef;
418
419	G4double epf1_in = 0.0;
420	G4double epf1_out = 0.0;
421
422	// C ----Matrice de rotation (noyau de fission -> Remnant)
423	if(peva >= 1.0e-4) {
424	sitet = std::sqrt(std::pow(pxeva,2)+std::pow(pyeva,2))/peva;
425	}
426	if(sitet > 1.0e-5) { //then
427	G4double cstet = pleva/peva;
428	G4double siphi = pyeva/(sitet*peva);
429	G4double csphi = pxeva/(sitet*peva);
430
431	R[1][1] = cstet*csphi;
432	R[1][2] = -siphi;
433	R[1][3] = sitet*csphi;
434	R[2][1] = cstet*siphi;
435	R[2][2] = csphi;
436	R[2][3] = sitet*siphi;
437	R[3][1] = -sitet;
438	R[3][2] = 0.0;
439	R[3][3] = cstet;
440	}
441	else {
442	R[1][1] = 1.0;
443	R[1][2] = 0.0;
444	R[1][3] = 0.0;
445	R[2][1] = 0.0;
446	R[2][2] = 1.0;
447	R[2][3] = 0.0;
448	R[3][1] = 0.0;
449	R[3][2] = 0.0;
450	R[3][3] = 1.0;
451	} // endif
452	// c test de verif:
453
454	if((zff1 <= 0.0) \|\| (aff1 <= 0.0) \|\| (aff1 < zff1)) { //then
455	if(verboseLevel > 2) {
456	G4cout <<"zf = " << zf <<" af = " << af <<"ee = " << ee <<"zff1 = " << zff1 <<"aff1 = " << aff1 << G4endl;
457	}
458	}
459	else {
460	// C ---------------------- PF1 will evaporate
461	epf1_in = double(eff1);
462	epf1_out = epf1_in;
463	// void evapora(G4double zprf, G4double aprf, G4double ee, G4double jprf,
464	// G4double zf_par, G4double af_par, G4double *mtota_par,
465	// G4double pleva_par, G4double pxeva_par, G4double *pyeva_par,
466	// G4double ff_par, G4int inttype_par, G4int *inum_par);
467	G4double zf1, af1, malpha1, ffpleva1, ffpxeva1, ffpyeva1;
468	G4int ff1, ftype1;
469	evapora(zff1, aff1, epf1_out, 0.0, &zf1, &af1, &malpha1, &ffpleva1,
470	&ffpxeva1, &ffpyeva1, &ff1, &ftype1, &inum);
471	// C On ajoute le fragment:
472	volant->iv = volant->iv + 1;
473	volant->acv[volant->iv] = af1;
474	volant->zpcv[volant->iv] = zf1;
475	if(verboseLevel > 2) {
476	G4cout << __FILE__ << ":" << __LINE__ << " Added: zf1 = " << zf1 << " af1 = " << af1 << " at index " << volant->iv << G4endl;
477	volant->dump();
478	}
479	if(verboseLevel > 2) {
480	G4cout <<"Added fission fragment: a = " << volant->acv[volant->iv] << " z = " << volant->zpcv[volant->iv] << G4endl;
481	}
482	peva = std::sqrt(std::pow(ffpxeva1,2) + std::pow(ffpyeva1,2) + std::pow(ffpleva1,2));
483	volant->pcv[volant->iv] = peva;
484	if(peva > 0.001) { // then
485	volant->xcv[volant->iv] = ffpxeva1/peva;
486	volant->ycv[volant->iv] = ffpyeva1/peva;
487	volant->zcv[volant->iv] = ffpleva1/peva;
488	}
489	else {
490	volant->xcv[volant->iv] = 1.0;
491	volant->ycv[volant->iv] = 0.0;
492	volant->zcv[volant->iv] = 0.0;
493	} // end if
494
495	// C Pour VÃ©rif evapo de PF1 dans le systeme du Noyau Fissionant
496	G4double bil1_e = 0.0;
497	G4double bil1_px = 0.0;
498	G4double bil1_py=0.0;
499	G4double bil1_pz=0.0;
500	for(G4int iloc = nbpevap + 1; iloc <= volant->iv; iloc++) { //do iloc=nbpevap+1,iv
501	// for(G4int iloc = nbpevap + 1; iloc <= volant->iv + 1; iloc++) { //do iloc=nbpevap+1,iv
502	mglms(volant->acv[iloc], volant->zpcv[iloc],0,&el);
503	masse = volant->zpcv[iloc]fmp + (volant->acv[iloc] - volant->zpcv[iloc])fmn + el;
504	bil1_e = bil1_e + std::sqrt(std::pow(volant->pcv[iloc],2) + std::pow(masse,2));
505	bil1_px = bil1_px + volant->pcv[iloc]*(volant->xcv[iloc]);
506	bil1_py = bil1_py + volant->pcv[iloc]*(volant->ycv[iloc]);
507	bil1_pz = bil1_pz + volant->pcv[iloc]*(volant->zcv[iloc]);
508	} // enddo
509
510	// Calcul des cosinus directeurs de PF1 dans le Remnant et calcul
511	// des coefs pour la transformation de Lorentz Systeme PF --> Systeme Remnant
512	translabpf(masse1,t1,p1,ctet1,phi1,gamfis,etfis,R,&plab1,&gam1,&eta1,csdir1);
513
514	// calcul des impulsions des particules evaporees dans le systeme Remnant:
515	if(verboseLevel > 2)
516	G4cout <<"2nd Translab (pf1 evap): Adding indices from " << nbpevap+1 << " to " << volant->iv << G4endl;
517	nopart = varntp->ntrack - 1;
518	translab(gam1,eta1,csdir1,nopart,nbpevap+1);
519	if(verboseLevel > 2) {
520	G4cout <<"After translab call... varntp->ntrack = " << varntp->ntrack << G4endl;
521	}
522	memiv = nbpevap + 1; // memoires pour la future transformation
523	mempaw = nopart; // remnant->labo pour pf1 et pf2.
524	lmi_pf1 = nopart + nbpevap + 1; // indices min et max dans /var_ntp/
525	lma_pf1 = nopart + volant->iv; // des particules issues de pf1
526	nbpevap = volant->iv; // nombre de particules d'evaporation traitees
527	} // end if
528	// C --------------------- End of PF1 calculation
529
530	// c test de verif:
531	if((zff2 <= 0.0) \|\| (aff2 <= 0.0) \|\| (aff2 <= zff2)) { //then
532	if(verboseLevel > 2) {
533	G4cout << zf << " " << af << " " << ee << " " << zff2 << " " << aff2 << G4endl;
534	}
535	}
536	else {
537	// C ---------------------- PF2 will evaporate
538	G4double epf2_in = double(eff2);
539	G4double epf2_out = epf2_in;
540	// void evapora(G4double zprf, G4double aprf, G4double ee, G4double jprf,
541	// G4double zf_par, G4double af_par, G4double *mtota_par,
542	// G4double pleva_par, G4double pxeva_par, G4double *pyeva_par,
543	// G4double ff_par, G4int inttype_par, G4int *inum_par);
544	G4double zf2, af2, malpha2, ffpleva2, ffpxeva2, ffpyeva2;
545	G4int ff2, ftype2;
546	evapora(zff2,aff2,epf2_out,0.0,&zf2,&af2,&malpha2,&ffpleva2,
547	&ffpxeva2,&ffpyeva2,&ff2,&ftype2,&inum);
548	// C On ajoute le fragment:
549	volant->iv = volant->iv + 1;
550	volant->acv[volant->iv] = af2;
551	volant->zpcv[volant->iv] = zf2;
552	if(verboseLevel > 2)
553	G4cout << __FILE__ << ":" << __LINE__ << " Added: zf2 = " << zf2 << " af2 = " << af2 << " at index " << volant->iv << G4endl;
554	if(verboseLevel > 2) {
555	G4cout <<"Added fission fragment: a = " << volant->acv[volant->iv] << " z = " << volant->zpcv[volant->iv] << G4endl;
556	}
557	peva = std::sqrt(std::pow(ffpxeva2,2) + std::pow(ffpyeva2,2) + std::pow(ffpleva2,2));
558	volant->pcv[volant->iv] = peva;
559	// exit(0);
560	if(peva > 0.001) { //then
561	volant->xcv[volant->iv] = ffpxeva2/peva;
562	volant->ycv[volant->iv] = ffpyeva2/peva;
563	volant->zcv[volant->iv] = ffpleva2/peva;
564	}
565	else {
566	volant->xcv[volant->iv] = 1.0;
567	volant->ycv[volant->iv] = 0.0;
568	volant->zcv[volant->iv] = 0.0;
569	} //end if
570	// C Pour VÃ©rif evapo de PF1 dans le systeme du Noyau Fissionant
571	G4double bil2_e = 0.0;
572	G4double bil2_px = 0.0;
573	G4double bil2_py = 0.0;
574	G4double bil2_pz = 0.0;
575	// for(G4int iloc = nbpevap + 1; iloc <= volant->iv; iloc++) { //do iloc=nbpevap+1,iv
576	for(G4int iloc = nbpevap + 1; iloc <= volant->iv; iloc++) { //do iloc=nbpevap+1,iv
577	mglms(volant->acv[iloc],volant->zpcv[iloc],0,&el);
578	masse = volant->zpcv[iloc]fmp + (volant->acv[iloc] - volant->zpcv[iloc])fmn + el;
579	bil2_e = bil2_e + std::sqrt(std::pow(volant->pcv[iloc],2) + std::pow(masse,2));
580	bil2_px = bil2_px + volant->pcv[iloc]*(volant->xcv[iloc]);
581	bil2_py = bil2_py + volant->pcv[iloc]*(volant->ycv[iloc]);
582	bil2_pz = bil2_pz + volant->pcv[iloc]*(volant->zcv[iloc]);
583	} //enddo
584
585	// C ----Calcul des cosinus directeurs de PF2 dans le Remnant et calcul
586	// c des coefs pour la transformation de Lorentz Systeme PF --> Systeme Remnant
587	G4double t2 = b - t1;
588	// G4double ctet2 = -ctet1;
589	ctet2 = -1.0*ctet1;
590	phi2 = dmod(phi1+3.141592654,6.283185308);
591	G4double p2 = std::sqrt(t2(t2+2.0masse2));
592
593	// void translabpf(G4double masse1, G4double t1, G4double p1, G4double ctet1,
594	// G4double phi1, G4double gamrem, G4double etrem, G4double R[][4],
595	// G4double plab1, G4double gam1, G4double *eta1, G4double csdir[]);
596	translabpf(masse2,t2,p2,ctet2,phi2,gamfis,etfis,R,&plab2,&gam2,&eta2,csdir2);
597	// C
598	// C calcul des impulsions des particules evaporees dans le systeme Remnant:
599	// c
600	if(verboseLevel > 2)
601	G4cout <<"3rd Translab (pf2 evap): Adding indices from " << nbpevap+1 << " to " << volant->iv << G4endl;
602	nopart = varntp->ntrack - 1;
603	translab(gam2,eta2,csdir2,nopart,nbpevap+1);
604	lmi_pf2 = nopart + nbpevap + 1; // indices min et max dans /var_ntp/
605	lma_pf2 = nopart + volant->iv; // des particules issues de pf2
606	} // end if
607	// C --------------------- End of PF2 calculation
608
609	// C Pour vÃ©rifications: calculs du noyau fissionant et des PF dans
610	// C le systeme du remnant.
611	for(G4int iloc = 1; iloc <= 3; iloc++) { // do iloc=1,3
612	pfis_rem[iloc] = 0.0;
613	} // enddo
614	G4double efis_rem, pfis_trav[4];
615	lorab(gamfis,etfis,massef,pfis_rem,&efis_rem,pfis_trav);
616	rotab(R,pfis_trav,pfis_rem);
617
618	stet1 = std::sqrt(1.0 - std::pow(ctet1,2));
619	pf1_rem[1] = p1stet1std::cos(phi1);
620	pf1_rem[2] = p1stet1std::sin(phi1);
621	pf1_rem[3] = p1*ctet1;
622	G4double e1_rem;
623	lorab(gamfis,etfis,masse1+t1,pf1_rem,&e1_rem,pfis_trav);
624	rotab(R,pfis_trav,pf1_rem);
625
626	stet2 = std::sqrt(1.0 - std::pow(ctet2,2));
627
628	G4double pf2_rem[4];
629	G4double e2_rem;
630	pf2_rem[1] = p2stet2std::cos(phi2);
631	pf2_rem[2] = p2stet2std::sin(phi2);
632	pf2_rem[3] = p2*ctet2;
633	lorab(gamfis,etfis,masse2+t2,pf2_rem,&e2_rem,pfis_trav);
634	rotab(R,pfis_trav,pf2_rem);
635	// C Verif 0: Remnant = evapo_pre_fission + Noyau Fissionant
636	bil_e = remmass - efis_rem - bil_e;
637	bil_px = bil_px + pfis_rem[1];
638	bil_py = bil_py + pfis_rem[2];
639	bil_pz = bil_pz + pfis_rem[3];
640	// C Verif 1: noyau fissionant = PF1 + PF2 dans le systeme remnant
641	// G4double bilan_e = efis_rem - e1_rem - e2_rem;
642	// G4double bilan_px = pfis_rem[1] - pf1_rem[1] - pf2_rem[1];
643	// G4double bilan_py = pfis_rem[2] - pf1_rem[2] - pf2_rem[2];
644	// G4double bilan_pz = pfis_rem[3] - pf1_rem[3] - pf2_rem[3];
645	// C Verif 2: PF1 et PF2 egaux a toutes leurs particules evaporees
646	// C (Systeme remnant)
647	if((lma_pf1-lmi_pf1) != 0) { //then
648	G4double bil_e_pf1 = e1_rem - epf1_out;
649	G4double bil_px_pf1 = pf1_rem[1];
650	G4double bil_py_pf1 = pf1_rem[2];
651	G4double bil_pz_pf1 = pf1_rem[3];
652	for(G4int ipf1 = lmi_pf1; ipf1 <= lma_pf1; ipf1++) { //do ipf1=lmi_pf1,lma_pf1
653	if(varntp->enerj[ipf1] <= 0.0) continue; // Safeguard against a division by zero
654	bil_e_pf1 = bil_e_pf1 - (std::pow(varntp->plab[ipf1],2) + std::pow(varntp->enerj[ipf1],2))/(2.0*(varntp->enerj[ipf1]));
655	cst = std::cos(varntp->tetlab[ipf1]/57.2957795);
656	sst = std::sin(varntp->tetlab[ipf1]/57.2957795);
657	csf = std::cos(varntp->philab[ipf1]/57.2957795);
658	ssf = std::sin(varntp->philab[ipf1]/57.2957795);
659	bil_px_pf1 = bil_px_pf1 - varntp->plab[ipf1]sstcsf;
660	bil_py_pf1 = bil_py_pf1 - varntp->plab[ipf1]sstssf;
661	bil_pz_pf1 = bil_pz_pf1 - varntp->plab[ipf1]*cst;
662	} // enddo
663	} //endif
664
665	if((lma_pf2-lmi_pf2) != 0) { //then
666	G4double bil_e_pf2 = e2_rem - epf2_out;
667	G4double bil_px_pf2 = pf2_rem[1];
668	G4double bil_py_pf2 = pf2_rem[2];
669	G4double bil_pz_pf2 = pf2_rem[3];
670	for(G4int ipf2 = lmi_pf2; ipf2 <= lma_pf2; ipf2++) { //do ipf2=lmi_pf2,lma_pf2
671	if(varntp->enerj[ipf2] <= 0.0) continue; // Safeguard against a division by zero
672	bil_e_pf2 = bil_e_pf2 - (std::pow(varntp->plab[ipf2],2) + std::pow(varntp->enerj[ipf2],2))/(2.0*(varntp->enerj[ipf2]));
673	G4double cst = std::cos(varntp->tetlab[ipf2]/57.2957795);
674	G4double sst = std::sin(varntp->tetlab[ipf2]/57.2957795);
675	G4double csf = std::cos(varntp->philab[ipf2]/57.2957795);
676	G4double ssf = std::sin(varntp->philab[ipf2]/57.2957795);
677	bil_px_pf2 = bil_px_pf2 - varntp->plab[ipf2]sstcsf;
678	bil_py_pf2 = bil_py_pf2 - varntp->plab[ipf2]sstssf;
679	bil_pz_pf2 = bil_pz_pf2 - varntp->plab[ipf2]*cst;
680	} // enddo
681	} //endif
682	// C
683	// C ---- Transformation systeme Remnant -> systeme labo. (evapo de PF1 ET PF2)
684	// C
685	// G4double mempaw, memiv;
686	if(verboseLevel > 2)
687	G4cout <<"4th Translab: Adding indices from " << memiv << " to " << volant->iv << G4endl;
688	translab(gamrem,etrem,csrem,mempaw,memiv);
689	// C ***************** END of fission calculations **********************
690	if(verboseLevel > 2) {
691	G4cout <<"Dump at the end of fission event " << G4endl;
692	volant->dump();
693	G4cout <<"End of dump." << G4endl;
694	}
695	}
696	else {
697	// C ********************** Evapo sans fission ***************************
698	// C Here, FF=0, --> Evapo sans fission, on ajoute le fragment:
699	// C *************************************************************************
700	varntp->kfis = 0;
701	if(verboseLevel > 2) {
702	G4cout <<"Evaporation without fission" << G4endl;
703	}
704	volant->iv = volant->iv + 1;
705	volant->acv[volant->iv] = af;
706	volant->zpcv[volant->iv] = zf;
707	if(verboseLevel > 2)
708	G4cout << __FILE__ << ":" << __LINE__ << " Added: zf = " << zf << " af = " << af << " at index " << volant->iv << G4endl;
709	G4double peva = std::sqrt(std::pow(pxeva,2)+std::pow(pyeva,2)+std::pow(pleva,2));
710	volant->pcv[volant->iv] = peva;
711	if(peva > 0.001) { //then
712	volant->xcv[volant->iv] = pxeva/peva;
713	volant->ycv[volant->iv] = pyeva/peva;
714	volant->zcv[volant->iv] = pleva/peva;
715	}
716	else {
717	volant->xcv[volant->iv] = 1.0;
718	volant->ycv[volant->iv] = 0.0;
719	volant->zcv[volant->iv] = 0.0;
720	} // end if
721
722	// C
723	// C calcul des impulsions des particules evaporees dans le systeme labo:
724	// c
725	trem = double(erecrem);
726	// C REMMASS = pace2(APRF,ZPRF) + APRFUMA - ZPRFMELEC !Canonic
727	// C REMMASS = MCOREM + DBLE(ESREM) !OK
728	remmass = mcorem; //Cugnon
729	// C GAMREM = (REMMASS + TREM)/REMMASS !OK
730	// C ETREM = DSQRT(TREM(TREM + 2.REMMASS))/REMMASS !OK
731	csrem[0] = 0.0; // Should not be used.
732	csrem[1] = alrem;
733	csrem[2] = berem;
734	csrem[3] = garem;
735
736	// for(G4int j = 1; j <= volant->iv; j++) { //do j=1,iv
737	for(G4int j = 1; j <= volant->iv; j++) { //do j=1,iv
738	if(volant->acv[j] == 0) {
739	if(verboseLevel > 2) {
740	G4cout <<"volant->acv[" << j << "] = 0" << G4endl;
741	G4cout <<"volant->iv = " << volant->iv << G4endl;
742	}
743	}
744	if(volant->acv[j] > 0) {
745	mglms(volant->acv[j],volant->zpcv[j],0,&el);
746	fmcv = volant->zpcv[j]fmp + (volant->acv[j] - volant->zpcv[j])fmn + el;
747	e_evapo = e_evapo + std::sqrt(std::pow(volant->pcv[j],2) + std::pow(fmcv,2));
748	}
749	} // enddo
750
751	// C Redefinition pour conservation d'impulsion!!!
752	// C this mass obtained by energy balance is very close to the
753	// C mass of the remnant computed by pace2 + excitation energy (EE). (OK)
754	remmass = e_evapo;
755
756	G4double gamrem = std::sqrt(std::pow(pcorem,2)+std::pow(remmass,2))/remmass;
757	G4double etrem = pcorem/remmass;
758
759	if(verboseLevel > 2)
760	G4cout <<"5th Translab (no fission): Adding indices from " << 1 << " to " << volant->iv << G4endl;
761	nopart = varntp->ntrack - 1;
762	translab(gamrem,etrem,csrem,nopart,1);
763
764	if(verboseLevel > 2) {
765	G4cout <<"Dump at the end of evaporation event " << G4endl;
766	volant->dump();
767	G4cout <<"End of dump." << G4endl;
768	}
769	// C End of the (FISSION - NO FISSION) condition (FF=1 or 0)
770	} //end if
771	// C ********************* FIN de l'EVAPO KHS ******************
772	}
773
774	// Evaporation code
775	void G4Abla::initEvapora()
776	{
777	// 37 C PROJECTILE AND TARGET PARAMETERS + CROSS SECTIONS
778	// 38 C COMMON /ABLAMAIN/ AP,ZP,AT,ZT,EAP,BETA,BMAXNUC,CRTOT,CRNUC,
779	// 39 C R_0,R_P,R_T, IMAX,IRNDM,PI,
780	// 40 C BFPRO,SNPRO,SPPRO,SHELL
781	// 41 C
782	// 42 C AP,ZP,AT,ZT - PROJECTILE AND TARGET MASSES
783	// 43 C EAP,BETA - BEAM ENERGY PER NUCLEON, V/C
784	// 44 C BMAXNUC - MAX. IMPACT PARAMETER FOR NUCL. REAC.
785	// 45 C CRTOT,CRNUC - TOTAL AND NUCLEAR REACTION CROSS SECTION
786	// 46 C R_0,R_P,R_T, - RADIUS PARAMETER, PROJECTILE+ TARGET RADII
787	// 47 C IMAX,IRNDM,PI - MAXIMUM NUMBER OF EVENTS, DUMMY, 3.141...
788	// 48 C BFPRO - FISSION BARRIER OF THE PROJECTILE
789	// 49 C SNPRO - NEUTRON SEPARATION ENERGY OF THE PROJECTILE
790	// 50 C SPPRO - PROTON " " " " "
791	// 51 C SHELL - GROUND STATE SHELL CORRECTION
792	// 52 C---------------------------------------------------------------------
793	// 53 C
794	// 54 C ENERGIES WIDTHS AND CROSS SECTIONS FOR EM EXCITATION
795	// 55 C COMMON /EMDPAR/ EGDR,EGQR,FWHMGDR,FWHMGQR,CREMDE1,CREMDE2,
796	// 56 C AE1,BE1,CE1,AE2,BE2,CE2,SR1,SR2,XR
797	// 57 C
798	// 58 C EGDR,EGQR - MEAN ENERGY OF GDR AND GQR
799	// 59 C FWHMGDR,FWHMGQR - FWHM OF GDR, GQR
800	// 60 C CREMDE1,CREMDE2 - EM CROSS SECTION FOR E1 AND E2
801	// 61 C AE1,BE1,CE1 - ARRAYS TO CALCULATE
802	// 62 C AE2,BE2,CE2 - THE EXCITATION ENERGY AFTER E.M. EXC.
803	// 63 C SR1,SR2,XR - WITH MONTE CARLO
804	// 64 C---------------------------------------------------------------------
805	// 65 C
806	// 66 C DEFORMATIONS AND G.S. SHELL EFFECTS
807	// 67 C COMMON /ECLD/ ECGNZ,ECFNZ,VGSLD,ALPHA
808	// 68 C
809	// 69 C ECGNZ - GROUND STATE SHELL CORR. FRLDM FOR A SPHERICAL G.S.
810	// 70 C ECFNZ - SHELL CORRECTION FOR THE SADDLE POINT (NOW: == 0)
811	// 71 C VGSLD - DIFFERENCE BETWEEN DEFORMED G.S. AND LDM VALUE
812	// 72 C ALPHA - ALPHA GROUND STATE DEFORMATION (THIS IS NOT BETA2!)
813	// 73 C BETA2 = SQRT(5/(4PI)) * ALPHA
814	// 74 C---------------------------------------------------------------------
815	// 75 C
816	// 76 C ARRAYS FOR EXCITATION ENERGY BY STATISTICAL HOLE ENERY MODEL
817	// 77 C COMMON /EENUC/ SHE, XHE
818	// 78 C
819	// 79 C SHE, XHE - ARRAYS TO CALCULATE THE EXC. ENERGY AFTER
820	// 80 C ABRASION BY THE STATISTICAL HOLE ENERGY MODEL
821	// 81 C---------------------------------------------------------------------
822	// 82 C
823	// 83 C G.S. SHELL EFFECT
824	// 84 C COMMON /EC2SUB/ ECNZ
825	// 85 C
826	// 86 C ECNZ G.S. SHELL EFFECT FOR THE MASSES (IDENTICAL TO ECGNZ)
827	// 87 C---------------------------------------------------------------------
828	// 88 C
829	// 89 C OPTIONS AND PARAMETERS FOR FISSION CHANNEL
830	// 90 C COMMON /FISS/ AKAP,BET,HOMEGA,KOEFF,IFIS,
831	// 91 C OPTSHP,OPTXFIS,OPTLES,OPTCOL
832	// 92 C
833	// 93 C AKAP - HBAR*2/(2 MN * R_0**2) = 10 MEV
834	// 94 C BET - REDUCED NUCLEAR FRICTION COEFFICIENT IN (1021 S-1)
835	// 95 C HOMEGA - CURVATURE OF THE FISSION BARRIER = 1 MEV
836	// 96 C KOEFF - COEFFICIENT FOR THE LD FISSION BARRIER == 1.0
837	// 97 C IFIS - 0/1 FISSION CHANNEL OFF/ON
838	// 98 C OPTSHP - INTEGER SWITCH FOR SHELL CORRECTION IN MASSES/ENERGY
839	// 99 C = 0 NO MICROSCOPIC CORRECTIONS IN MASSES AND ENERGY
840	// 100 C = 1 SHELL , NO PAIRING
841	// 101 C = 2 PAIRING, NO SHELL
842	// 102 C = 3 SHELL AND PAIRING
843	// 103 C OPTCOL - 0/1 COLLECTIVE ENHANCEMENT SWITCHED ON/OFF
844	// 104 C OPTXFIS- 0,1,2 FOR MYERS & SWIATECKI, DAHLINGER, ANDREYEV
845	// 105 C FISSILITY PARAMETER.
846	// 106 C OPTLES - CONSTANT TEMPERATURE LEVEL DENSITY FOR A,Z > TH-224
847	// 107 C OPTCOL - 0/1 COLLECTIVE ENHANCEMENT OFF/ON
848	// 108 C---------------------------------------------------------------------
849	// 109 C
850	// 110 C OPTIONS
851	// 111 C COMMON /OPT/ OPTEMD,OPTCHA,EEFAC
852	// 112 C
853	// 113 C OPTEMD - 0/1 NO EMD / INCL. EMD
854	// 114 C OPTCHA - 0/1 0 GDR / 1 HYPERGEOMETRICAL PREFRAGMENT-CHARGE-DIST.
855	// 115 C * RECOMMENDED IS OPTCHA = 1 *
856	// 116 C EEFAC - EXCITATION ENERGY FACTOR, 2.0 RECOMMENDED
857	// 117 C---------------------------------------------------------------------
858	// 118 C
859	// 119 C FISSION BARRIERS
860	// 120 C COMMON /FB/ EFA
861	// 121 C EFA - ARRAY OF FISSION BARRIERS
862	// 122 C---------------------------------------------------------------------
863	// 123 C
864	// 124 C p LEVEL DENSITY PARAMETERS
865	// 125 C COMMON /ALD/ AV,AS,AK,OPTAFAN
866	// 126 C AV,AS,AK - VOLUME,SURFACE,CURVATURE DEPENDENCE OF THE
867	// 127 C LEVEL DENSITY PARAMETER
868	// 128 C OPTAFAN - 0/1 AF/AN >=1 OR AF/AN ==1
869	// 129 C RECOMMENDED IS OPTAFAN = 0
870	// 130 C---------------------------------------------------------------------
871	// 131 C ____________________________________________________________________
872	// 132 C /
873	// 133 C / INITIALIZES PARAMETERS IN COMMON /ABRAMAIN/, /EMDPAR/, /ECLD/ ...
874	// 134 C / PROJECTILE PARAMETERS, EMD PARAMETERS, SHELL CORRECTION TABLES.
875	// 135 C / CALCULATES MAXIMUM IMPACT PARAMETER FOR NUCLEAR COLLISIONS AND
876	// 136 C / TOTAL GEOMETRICAL CROSS SECTION + EMD CROSS SECTIONS
877	// 137 C ____________________________________________________________________
878	// 138 C
879	// 139 C
880	// 201 C
881	// 202 C---------- SET INPUT VALUES
882	// 203 C
883	// 204 C *** INPUT FROM UNIT 10 IN THE FOLLOWING SEQUENCE !
884	// 205 C AP1 = INTEGER !
885	// 206 C ZP1 = INTEGER !
886	// 207 C AT1 = INTEGER !
887	// 208 C ZT1 = INTEGER !
888	// 209 C EAP = REAL !
889	// 210 C IMAX = INTEGER !
890	// 211 C IFIS = INTEGER SWITCH FOR FISSION
891	// 212 C OPTSHP = INTEGER SWITCH FOR SHELL CORRECTION IN MASSES/ENERGY
892	// 213 C =0 NO MICROSCOPIC CORRECTIONS IN MASSES AND ENERGY
893	// 214 C =1 SHELL , NO PAIRING CORRECTION
894	// 215 C =2 PAIRING, NO SHELL CORRECTION
895	// 216 C =3 SHELL AND PAIRING CORRECTION IN MASSES AND ENERGY
896	// 217 C OPTEMD =0,1 0 NO EMD, 1 INCL. EMD
897	// 218 C ELECTROMAGNETIC DISSOZIATION IS CALCULATED AS WELL.
898	// 219 C OPTCHA =0,1 0 GDR- , 1 HYPERGEOMETRICAL PREFRAGMENT-CHARGE-DIST.
899	// 220 C RECOMMENDED IS OPTCHA=1
900	// 221 C OPTCOL =0,1 COLLECTIVE ENHANCEMENT SWITCHED ON 1 OR OFF 0 IN DENSN
901	// 222 C OPTAFAN=0,1 SWITCH FOR AF/AN = 1 IN DENSNIV 0 AF/AN>1 1 AF/AN=1
902	// 223 C AKAP = REAL ALWAYS EQUALS 10
903	// 224 C BET = REAL REDUCED FRICTION COEFFICIENT / 10(+21) S(-1)
904	// 225 C HOMEGA = REAL CURVATURE / MEV RECOMMENDED = 1. MEV
905	// 226 C KOEFF = REAL COEFFICIENT FOR FISSION BARRIER
906	// 227 C OPTXFIS= INTEGER 0,1,2 FOR MYERS & SWIATECKI, DAHLINGER, ANDREYEV
907	// 228 C FISSILITY PARAMETER.
908	// 229 C EEFAC = REAL EMPIRICAL FACTOR FOR THE EXCITATION ENERGY
909	// 230 C RECOMMENDED 2.D0, STATISTICAL ABRASION MODELL 1.D0
910	// 231 C AV = REAL KOEFFICIENTS FOR CALCULATION OF A(TILDE)
911	// 232 C AS = REAL LEVEL DENSITY PARAMETER
912	// 233 C AK = REAL
913	// 234 C
914	// 235 C This following inputs will be initialized in the main through the
915	// 236 C common /ABLAMAIN/ (A.B.)
916	// 237
917
918	// switch-fission.1=on.0=off
919	fiss->ifis = 1;
920
921	// shell+pairing.0-1-2-3
922	fiss->optshp = 0;
923
924	// optemd =0,1 0 no emd, 1 incl. emd
925	opt->optemd = 1;
926	// read(10,*,iostat=io) dum(10),optcha
927	opt->optcha = 1;
928
929	// not.to.be.changed.(akap)
930	fiss->akap = 10.0;
931
932	// nuclear.viscosity.(beta)
933	fiss->bet = 1.5;
934
935	// potential-curvature
936	fiss->homega = 1.0;
937
938	// fission-barrier-coefficient
939	fiss->koeff = 1.;
940
941	//collective enhancement switched on 1 or off 0 in densn (qr=val or =1.)
942	fiss->optcol = 0;
943
944	// switch-for-low-energy-sys
945	fiss->optles = 0;
946
947	opt->eefac = 2.;
948
949	ald->optafan = 0;
950
951	ald->av = 0.073e0;
952	ald->as = 0.095e0;
953	ald->ak = 0.0e0;
954
955	if(verboseLevel > 3) {
956	G4cout <<"ifis " << fiss->ifis << G4endl;
957	G4cout <<"optshp " << fiss->optshp << G4endl;
958	G4cout <<"optemd " << opt->optemd << G4endl;
959	G4cout <<"optcha " << opt->optcha << G4endl;
960	G4cout <<"akap " << fiss->akap << G4endl;
961	G4cout <<"bet " << fiss->bet << G4endl;
962	G4cout <<"homega " << fiss->homega << G4endl;
963	G4cout <<"koeff " << fiss->koeff << G4endl;
964	G4cout <<"optcol " << fiss->optcol << G4endl;
965	G4cout <<"optles " << fiss->optles << G4endl;
966	G4cout <<"eefac " << opt->eefac << G4endl;
967	G4cout <<"optafan " << ald->optafan << G4endl;
968	G4cout <<"av " << ald->av << G4endl;
969	G4cout <<"as " << ald->as << G4endl;
970	G4cout <<"ak " << ald->ak << G4endl;
971	}
972	fiss->optxfis = 1;
973
974	G4InclAblaDataFile *dataInterface = new G4InclAblaDataFile();
975	if(dataInterface->readData() == true) {
976	if(verboseLevel > 0) {
977	G4cout <<"G4Abla: Datafiles read successfully." << G4endl;
978	}
979	}
980	else {
981	G4Exception("ERROR: Failed to read datafiles.");
982	}
983
984	for(int z = 0; z < 99; z++) { //do 30 z = 0,98,1
985	for(int n = 0; n < 154; n++) { //do 31 n = 0,153,1
986	ecld->ecfnz[n][z] = 0.e0;
987	ec2sub->ecnz[n][z] = dataInterface->getEcnz(n,z);
988	ecld->ecgnz[n][z] = dataInterface->getEcnz(n,z);
989	ecld->alpha[n][z] = dataInterface->getAlpha(n,z);
990	ecld->vgsld[n][z] = dataInterface->getVgsld(n,z);
991	// if(ecld->ecgnz[n][z] != 0.0) G4cout <<"ecgnz[" << n << "][" << z << "] = " << ecld->ecgnz[n][z] << G4endl;
992	}
993	}
994
995	for(int z = 0; z < 500; z++) {
996	for(int a = 0; a < 500; a++) {
997	pace->dm[z][a] = dataInterface->getPace2(z,a);
998	}
999	}
1000
1001	delete dataInterface;
1002	}
1003
1004	void G4Abla::qrot(G4double z, G4double a, G4double bet, G4double sig, G4double u, G4double *qr)
1005	{
1006	G4double ucr = 10.0; // Critical energy for damping.
1007	G4double dcr = 40.0; // Width of damping.
1008	G4double ponq = 0.0, dn = 0.0, n = 0.0, dz = 0.0;
1009
1010	if(((std::fabs(bet)-1.15) < 0) \|\| ((std::fabs(bet)-1.15) == 0)) {
1011	goto qrot10;
1012	}
1013
1014	if((std::fabs(bet)-1.15) > 0) {
1015	goto qrot11;
1016	}
1017
1018	qrot10:
1019	n = a - z;
1020	dz = std::fabs(z - 82.0);
1021	if (n > 104) {
1022	dn = std::fabs(n-126.e0);
1023	}
1024	else {
1025	dn = std::fabs(n - 82.0);
1026	}
1027
1028	bet = 0.022 + 0.003dn + 0.005dz;
1029
1030	sig = 25.0std::pow(bet,2) sig;
1031
1032	qrot11:
1033	ponq = (u - ucr)/dcr;
1034
1035	if (ponq > 700.0) {
1036	ponq = 700.0;
1037	}
1038	if (sig < 1.0) {
1039	sig = 1.0;
1040	}
1041	(qr) = 1.0/(1.0 + std::exp(ponq)) (sig - 1.0) + 1.0;
1042
1043	if ((*qr) < 1.0) {
1044	(*qr) = 1.0;
1045	}
1046
1047	return;
1048	}
1049
1050	void G4Abla::mglw(G4double a, G4double z, G4double *el)
1051	{
1052	// MODEL DE LA GOUTTE LIQUIDE DE C. F. WEIZSACKER.
1053	// USUALLY AN OBSOLETE OPTION
1054
1055	G4int a1 = 0, z1 = 0;
1056	G4double xv = 0.0, xs = 0.0, xc = 0.0, xa = 0.0;
1057
1058	a1 = idnint(a);
1059	z1 = idnint(z);
1060
1061	if ((a <= 0.01) \|\| (z < 0.01)) {
1062	(*el) = 1.0e38;
1063	}
1064	else {
1065	xv = -15.56*a;
1066	xs = 17.23*std::pow(a,(2.0/3.0));
1067
1068	if (a > 1.0) {
1069	xc = 0.7z(z-1.0)*std::pow((a-1.0),(-1.e0/3.e0));
1070	}
1071	else {
1072	xc = 0.0;
1073	}
1074	}
1075
1076	xa = 23.6(std::pow((a-2.0z),2)/a);
1077	(*el) = xv+xs+xc+xa;
1078	return;
1079	}
1080
1081	void G4Abla::mglms(G4double a, G4double z, G4int refopt4, G4double *el)
1082	{
1083	// USING FUNCTION EFLMAC(IA,IZ,0)
1084	//
1085	// REFOPT4 = 0 : WITHOUT MICROSCOPIC CORRECTIONS
1086	// REFOPT4 = 1 : WITH SHELL CORRECTION
1087	// REFOPT4 = 2 : WITH PAIRING CORRECTION
1088	// REFOPT4 = 3 : WITH SHELL- AND PAIRING CORRECTION
1089
1090	// 1839 C-----------------------------------------------------------------------
1091	// 1840 C A1 LOCAL MASS NUMBER (INTEGER VARIABLE OF A)
1092	// 1841 C Z1 LOCAL NUCLEAR CHARGE (INTEGER VARIABLE OF Z)
1093	// 1842 C REFOPT4 OPTION, SPECIFYING THE MASS FORMULA (SEE ABOVE)
1094	// 1843 C A MASS NUMBER
1095	// 1844 C Z NUCLEAR CHARGE
1096	// 1845 C DEL PAIRING CORRECTION
1097	// 1846 C EL BINDING ENERGY
1098	// 1847 C ECNZ( , ) TABLE OF SHELL CORRECTIONS
1099	// 1848 C-----------------------------------------------------------------------
1100	// 1849 C
1101	G4int a1 = idnint(a);
1102	G4int z1 = idnint(z);
1103
1104	if ( (a1 <= 0) \|\| (z1 <= 0) \|\| ((a1-z1) <= 0) ) { //then
1105	// modif pour rÃ©cupÃ©rer une masse p et n correcte:
1106	(*el) = 0.0;
1107	return;
1108	// goto mglms50;
1109	}
1110	else {
1111	// binding energy incl. pairing contr. is calculated from
1112	// function eflmac
1113	(*el) = eflmac(a1,z1,0,refopt4);
1114	if (refopt4 > 0) {
1115	if (refopt4 != 2) {
1116	(el) = (el) + ec2sub->ecnz[a1-z1][z1];
1117	//(el) = (el) + ec2sub->ecnz[z1][a1-z1];
1118	}
1119	}
1120	}
1121	return;
1122	}
1123
1124	G4double G4Abla::spdef(G4int a, G4int z, G4int optxfis)
1125	{
1126
1127	// INPUT: A,Z,OPTXFIS MASS AND CHARGE OF A NUCLEUS,
1128	// OPTION FOR FISSILITY
1129	// OUTPUT: SPDEF
1130
1131	// ALPHA2 SADDLE POINT DEF. COHEN&SWIATECKI ANN.PHYS. 22 (1963) 406
1132	// RANGING FROM FISSILITY X=0.30 TO X=1.00 IN STEPS OF 0.02
1133
1134	G4int index = 0;
1135	G4double x = 0.0, v = 0.0, dx = 0.0;
1136
1137	const G4int alpha2Size = 37;
1138	// The value 0.0 at alpha2[0] added by PK.
1139	G4double alpha2[alpha2Size] = {0.0, 2.5464e0, 2.4944e0, 2.4410e0, 2.3915e0, 2.3482e0,
1140	2.3014e0, 2.2479e0, 2.1982e0, 2.1432e0, 2.0807e0, 2.0142e0, 1.9419e0,
1141	1.8714e0, 1.8010e0, 1.7272e0, 1.6473e0, 1.5601e0, 1.4526e0, 1.3164e0,
1142	1.1391e0, 0.9662e0, 0.8295e0, 0.7231e0, 0.6360e0, 0.5615e0, 0.4953e0,
1143	0.4354e0, 0.3799e0, 0.3274e0, 0.2779e0, 0.2298e0, 0.1827e0, 0.1373e0,
1144	0.0901e0, 0.0430e0, 0.0000e0};
1145
1146	dx = 0.02;
1147	x = fissility(a,z,optxfis);
1148
1149	if (x > 1.0) {
1150	x = 1.0;
1151	}
1152
1153	if (x < 0.0) {
1154	x = 0.0;
1155	}
1156
1157	v = (x - 0.3)/dx + 1.0;
1158	index = idnint(v);
1159
1160	if (index < 1) {
1161	return(alpha2[1]); // alpha2[0] -> alpha2[1]
1162	}
1163
1164	if (index == 36) { //then // :::FIXME:: Possible off-by-one bug...
1165	return(alpha2[36]);
1166	}
1167	else {
1168	return(alpha2[index] + (alpha2[index+1] - alpha2[index]) / dx * ( x - (0.3e0 + dx*(index-1)))); //:::FIXME::: Possible off-by-one
1169	}
1170
1171	return alpha2[0]; // The algorithm is not supposed to reach this point.
1172	}
1173
1174	G4double G4Abla::fissility(int a,int z, int optxfis)
1175	{
1176	// CALCULATION OF FISSILITY PARAMETER
1177	//
1178	// INPUT: A,Z INTEGER MASS & CHARGE OF NUCLEUS
1179	// OPTXFIS = 0 : MYERS, SWIATECKI
1180	// 1 : DAHLINGER
1181	// 2 : ANDREYEV
1182
1183	G4double aa = 0.0, zz = 0.0, i = 0.0;
1184	G4double fissilityResult = 0.0;
1185
1186	aa = double(a);
1187	zz = double(z);
1188	i = double(a-2*z) / aa;
1189
1190	// myers & swiatecki droplet modell
1191	if (optxfis == 0) { //then
1192	fissilityResult = std::pow(zz,2) / aa /50.8830e0 / (1.0e0 - 1.7826e0 * std::pow(i,2));
1193	}
1194
1195	if (optxfis == 1) {
1196	// dahlinger fit:
1197	fissilityResult = std::pow(zz,2) / aa * std::pow((49.22e0(1.e0 - 0.3803e0std::pow(i,2) - 20.489e0*std::pow(i,4))),(-1));
1198	}
1199
1200	if (optxfis == 2) {
1201	// dubna fit:
1202	fissilityResult = std::pow(zz,2) / aa /(48.e0(1.e0 - 17.22e0std::pow(i,4)));
1203	}
1204
1205	return fissilityResult;
1206	}
1207
1208	void G4Abla::evapora(G4double zprf, G4double aprf, G4double ee, G4double jprf,
1209	G4double zf_par, G4double af_par, G4double *mtota_par,
1210	G4double pleva_par, G4double pxeva_par, G4double *pyeva_par,
1211	G4int ff_par, G4int inttype_par, G4int *inum_par)
1212	{
1213	G4double zf = (*zf_par);
1214	G4double af = (*af_par);
1215	G4double mtota = (*mtota_par);
1216	G4double pleva = (*pleva_par);
1217	G4double pxeva = (*pxeva_par);
1218	G4double pyeva = (*pyeva_par);
1219	G4int ff = (*ff_par);
1220	G4int inttype = (*inttype_par);
1221	G4int inum = (*inum_par);
1222
1223	G4int idebug = 0;
1224
1225	if(idebug == 666) {
1226	zprf = 81.;
1227	aprf = 201.;
1228	// ee = 86.5877686;
1229	ee = 300.0;
1230	jprf = 32.;
1231	zf = 0.;
1232	af = 0.;
1233	mtota = 0.;
1234	pleva = 0.;
1235	pxeva = 0.;
1236	pyeva = 0.;
1237	ff = -1;
1238	inttype = 0;
1239	inum = 2;
1240	G4cout <<" PK::: EVAPORA event " << inum << G4endl;
1241	G4cout <<" PK::: zprf = " << zprf << G4endl;
1242	G4cout <<" PK::: aprf = " << aprf << G4endl;
1243	G4cout <<" PK::: ee = " << ee << G4endl;
1244	G4cout <<" PK::: eeDiff = " << ee - 86.5877686 << G4endl;
1245	G4cout <<" PK::: jprf = " << jprf << G4endl;
1246	G4cout <<" PK::: zf = " << zf << G4endl;
1247	G4cout <<" PK::: af = " << af << G4endl;
1248	G4cout <<" PK::: mtota = " << mtota << G4endl;
1249	G4cout <<" PK::: pleva = " << pleva << G4endl;
1250	G4cout <<" PK::: pxeva = " << pxeva << G4endl;
1251	G4cout <<" PK::: pyeva = " << pyeva << G4endl;
1252	G4cout <<" PK::: ff = " << ff << G4endl;
1253	G4cout <<" PK::: inttype = " << inttype << G4endl;
1254	G4cout <<" PK::: inum = " << inum << G4endl;
1255	}
1256	// 533 C
1257	// 534 C INPUT:
1258	// 535 C
1259	// 536 C ZPRF, APRF, EE(EE IS MODIFIED!), JPRF
1260	// 537 C
1261	// 538 C PROJECTILE AND TARGET PARAMETERS + CROSS SECTIONS
1262	// 539 C COMMON /ABRAMAIN/ AP,ZP,AT,ZT,EAP,BETA,BMAXNUC,CRTOT,CRNUC,
1263	// 540 C R_0,R_P,R_T, IMAX,IRNDM,PI,
1264	// 541 C BFPRO,SNPRO,SPPRO,SHELL
1265	// 542 C
1266	// 543 C AP,ZP,AT,ZT - PROJECTILE AND TARGET MASSES
1267	// 544 C EAP,BETA - BEAM ENERGY PER NUCLEON, V/C
1268	// 545 C BMAXNUC - MAX. IMPACT PARAMETER FOR NUCL. REAC.
1269	// 546 C CRTOT,CRNUC - TOTAL AND NUCLEAR REACTION CROSS SECTION
1270	// 547 C R_0,R_P,R_T, - RADIUS PARAMETER, PROJECTILE+ TARGET RADII
1271	// 548 C IMAX,IRNDM,PI - MAXIMUM NUMBER OF EVENTS, DUMMY, 3.141...
1272	// 549 C BFPRO - FISSION BARRIER OF THE PROJECTILE
1273	// 550 C SNPRO - NEUTRON SEPARATION ENERGY OF THE PROJECTILE
1274	// 551 C SPPRO - PROTON " " " " "
1275	// 552 C SHELL - GROUND STATE SHELL CORRECTION
1276	// 553 C
1277	// 554 C---------------------------------------------------------------------
1278	// 555 C FISSION BARRIERS
1279	// 556 C COMMON /FB/ EFA
1280	// 557 C EFA - ARRAY OF FISSION BARRIERS
1281	// 558 C---------------------------------------------------------------------
1282	// 559 C OUTPUT:
1283	// 560 C ZF, AF, MTOTA, PLEVA, PTEVA, FF, INTTYPE, INUM
1284	// 561 C
1285	// 562 C ZF,AF - CHARGE AND MASS OF FINAL FRAGMENT AFTER EVAPORATION
1286	// 563 C MTOTA _ NUMBER OF EVAPORATED ALPHAS
1287	// 564 C PLEVA,PXEVA,PYEVA - MOMENTUM RECOIL BY EVAPORATION
1288	// 565 C INTTYPE - TYPE OF REACTION 0/1 NUCLEAR OR ELECTROMAGNETIC
1289	// 566 C FF - 0/1 NO FISSION / FISSION EVENT
1290	// 567 C INUM - EVENTNUMBER
1291	// 568 C ____________________________________________________________________
1292	// 569 C /
1293	// 570 C / CALCUL DE LA MASSE ET CHARGE FINALES D'UNE CHAINE D'EVAPORATION
1294	// 571 C /
1295	// 572 C / PROCEDURE FOR CALCULATING THE FINAL MASS AND CHARGE VALUES OF A
1296	// 573 C / SPECIFIC EVAPORATION CHAIN, STARTING POINT DEFINED BY (APRF, ZPRF,
1297	// 574 C / EE)
1298	// 575 C / On ajoute les 3 composantes de l'impulsion (PXEVA,PYEVA,PLEVA)
1299	// 576 C / (actuellement PTEVA n'est pas correct; mauvaise norme...)
1300	// 577 C /____________________________________________________________________
1301	// 578 C
1302	// 612 C
1303	// 613 C-----------------------------------------------------------------------
1304	// 614 C IRNDM DUMMY ARGUMENT FOR RANDOM-NUMBER FUNCTION
1305	// 615 C SORTIE LOCAL HELP VARIABLE TO END THE EVAPORATION CHAIN
1306	// 616 C ZF NUCLEAR CHARGE OF THE FRAGMENT
1307	// 617 C ZPRF NUCLEAR CHARGE OF THE PREFRAGMENT
1308	// 618 C AF MASS NUMBER OF THE FRAGMENT
1309	// 619 C APRF MASS NUMBER OF THE PREFRAGMENT
1310	// 620 C EPSILN ENERGY BURNED IN EACH EVAPORATION STEP
1311	// 621 C MALPHA LOCAL MASS CONTRIBUTION TO MTOTA IN EACH EVAPORATION
1312	// 622 C STEP
1313	// 623 C EE EXCITATION ENERGY (VARIABLE)
1314	// 624 C PROBP PROTON EMISSION PROBABILITY
1315	// 625 C PROBN NEUTRON EMISSION PROBABILITY
1316	// 626 C PROBA ALPHA-PARTICLE EMISSION PROBABILITY
1317	// 627 C PTOTL TOTAL EMISSION PROBABILITY
1318	// 628 C E LOWEST PARTICLE-THRESHOLD ENERGY
1319	// 629 C SN NEUTRON SEPARATION ENERGY
1320	// 630 C SBP PROTON SEPARATION ENERGY PLUS EFFECTIVE COULOMB
1321	// 631 C BARRIER
1322	// 632 C SBA ALPHA-PARTICLE SEPARATION ENERGY PLUS EFFECTIVE
1323	// 633 C COULOMB BARRIER
1324	// 634 C BP EFFECTIVE PROTON COULOMB BARRIER
1325	// 635 C BA EFFECTIVE ALPHA COULOMB BARRIER
1326	// 636 C MTOTA TOTAL MASS OF THE EVAPORATED ALPHA PARTICLES
1327	// 637 C X UNIFORM RANDOM NUMBER FOR NUCLEAR CHARGE
1328	// 638 C AMOINS LOCAL MASS NUMBER OF EVAPORATED PARTICLE
1329	// 639 C ZMOINS LOCAL NUCLEAR CHARGE OF EVAPORATED PARTICLE
1330	// 640 C ECP KINETIC ENERGY OF PROTON WITHOUT COULOMB
1331	// 641 C REPULSION
1332	// 642 C ECN KINETIC ENERGY OF NEUTRON
1333	// 643 C ECA KINETIC ENERGY OF ALPHA PARTICLE WITHOUT COULOMB
1334	// 644 C REPULSION
1335	// 645 C PLEVA TRANSVERSAL RECOIL MOMENTUM OF EVAPORATION
1336	// 646 C PTEVA LONGITUDINAL RECOIL MOMENTUM OF EVAPORATION
1337	// 647 C FF FISSION FLAG
1338	// 648 C INTTYPE INTERACTION TYPE FLAG
1339	// 649 C RNDX RECOIL MOMENTUM IN X-DIRECTION IN A SINGLE STEP
1340	// 650 C RNDY RECOIL MOMENTUM IN Y-DIRECTION IN A SINGLE STEP
1341	// 651 C RNDZ RECOIL MOMENTUM IN Z-DIRECTION IN A SINGLE STEP
1342	// 652 C RNDN NORMALIZATION OF RECOIL MOMENTUM FOR EACH STEP
1343	// 653 C-----------------------------------------------------------------------
1344	// 654 C
1345	// 655 SAVE
1346	// SAVE -> static
1347
1348	static G4int sortie = 0;
1349	static G4double epsiln = 0.0, probp = 0.0, probn = 0.0, proba = 0.0, ptotl = 0.0, e = 0.0;
1350	static G4double sn = 0.0, sbp = 0.0, sba = 0.0, x = 0.0, amoins = 0.0, zmoins = 0.0;
1351	G4double ecn = 0.0, ecp = 0.0,eca = 0.0, bp = 0.0, ba = 0.0;
1352	static G4double pteva = 0.0;
1353
1354	static G4int itest = 0;
1355	static G4double probf = 0.0;
1356
1357	static G4int k = 0, j = 0, il = 0;
1358
1359	static G4double ctet1 = 0.0, stet1 = 0.0, phi1 = 0.0;
1360	static G4double sbfis = 0.0, rnd = 0.0;
1361	static G4double selmax = 0.0;
1362	static G4double segs = 0.0;
1363	static G4double ef = 0.0;
1364	static G4int irndm = 0;
1365
1366	static G4double pc = 0.0, malpha = 0.0;
1367
1368	zf = zprf;
1369	af = aprf;
1370	pleva = 0.0;
1371	pteva = 0.0;
1372	pxeva = 0.0;
1373	pyeva = 0.0;
1374
1375	sortie = 0;
1376	ff = 0;
1377
1378	itest = 0;
1379	if (itest == 1) {
1380	G4cout << "***************************" << G4endl;
1381	}
1382
1383	evapora10:
1384
1385	if (itest == 1) {
1386	G4cout <<"------zf,af,ee------" << idnint(zf) << "," << idnint(af) << "," << ee << G4endl;
1387	}
1388
1389	// calculation of the probabilities for the different decay channels
1390	// plus separation energies and kinetic energies of the particles
1391	direct(zf,af,ee,jprf,&probp,&probn,&proba,&probf,&ptotl,
1392	&sn,&sbp,&sba,&ecn,&ecp,&eca,&bp,&ba,inttype,inum,itest); //:::FIXME::: Call
1393	if((eca+ba) < 0) {
1394	eca = 0.0;
1395	ba = 0.0;
1396	}
1397	k = idnint(zf);
1398	j = idnint(af-zf);
1399
1400	// now ef is calculated from efa that depends on the subroutine
1401	// barfit which takes into account the modification on the ang. mom.
1402	// jb mvr 6-aug-1999
1403	// note *** shell correction! (ecgnz) jb mvr 20-7-1999
1404	il = idnint(jprf);
1405	barfit(k,k+j,il,&sbfis,&segs,&selmax);
1406
1407	if ((fiss->optshp == 1) \|\| (fiss->optshp == 3)) { //then
1408	// fb->efa[k][j] = G4double(sbfis) - ecld->ecgnz[j][k];
1409	fb->efa[j][k] = G4double(sbfis) - ecld->ecgnz[j][k];
1410	}
1411	else {
1412	fb->efa[j][k] = G4double(sbfis);
1413	// fb->efa[j][k] = G4double(sbfis);
1414	} //end if
1415	ef = fb->efa[j][k];
1416	// ef = fb->efa[j][k];
1417	// here the final steps of the evaporation are calculated
1418	if ((sortie == 1) \|\| (ptotl == 0.e0)) {
1419	e = dmin1(sn,sbp,sba);
1420	if (e > 1.0e30) {
1421	if(verboseLevel > 2) {
1422	G4cout <<"erreur a la sortie evapora,e>1.e30,af=" << af <<" zf=" << zf << G4endl;
1423	}
1424	}
1425	if (zf <= 6.0) {
1426	goto evapora100;
1427	}
1428	if (e < 0.0) {
1429	if (sn == e) {
1430	af = af - 1.e0;
1431	}
1432	else if (sbp == e) {
1433	af = af - 1.0;
1434	zf = zf - 1.0;
1435	}
1436	else if (sba == e) {
1437	af = af - 4.0;
1438	zf = zf - 2.0;
1439	}
1440	if (af < 2.5) {
1441	goto evapora100;
1442	}
1443	goto evapora10;
1444	}
1445	goto evapora100;
1446	}
1447	irndm = irndm + 1;
1448
1449	// here the normal evaporation cascade starts
1450
1451	// random number for the evaporation
1452	// x = double(Rndm(irndm))*ptotl;
1453	// x = double(haz(1))*ptotl;
1454	x = randomGenerator->getRandom() * ptotl;
1455	// G4cout <<"proba = " << proba << G4endl;
1456	// G4cout <<"probp = " << probp << G4endl;
1457	// G4cout <<"probn = " << probn << G4endl;
1458	// G4cout <<"probf = " << probf << G4endl;
1459
1460	itest = 0;
1461	if (x < proba) {
1462	// alpha evaporation
1463	if (itest == 1) {
1464	G4cout <<"PK::: < alpha evaporation >" << G4endl;
1465	}
1466	amoins = 4.0;
1467	zmoins = 2.0;
1468	epsiln = sba + eca;
1469	pc = std::sqrt(std::pow((1.0 + (eca+ba)/3.72834e3),2) - 1.0) * 3.72834e3;
1470	malpha = 4.0;
1471
1472	// volant:
1473	volant->iv = volant->iv + 1;
1474	volant->acv[volant->iv] = 4.;
1475	volant->zpcv[volant->iv] = 2.;
1476	volant->pcv[volant->iv] = pc;
1477	}
1478	else if (x < proba+probp) {
1479	// proton evaporation
1480	if (itest == 1) {
1481	G4cout <<"PK::: < proton evaporation >" << G4endl;
1482	}
1483	amoins = 1.0;
1484	zmoins = 1.0;
1485	epsiln = sbp + ecp;
1486	pc = std::sqrt(std::pow((1.0 + (ecp + bp)/9.3827e2),2) - 1.0) * 9.3827e2;
1487	malpha = 0.0;
1488	// volant:
1489	volant->iv = volant->iv + 1;
1490	volant->acv[volant->iv] = 1.0;
1491	volant->zpcv[volant->iv] = 1.;
1492	volant->pcv[volant->iv] = pc;
1493	}
1494	else if (x < proba+probp+probn) {
1495	// neutron evaporation
1496	if (itest == 1) {
1497	G4cout <<"PK::: < neutron evaporation >" << G4endl;
1498	}
1499	amoins = 1.0;
1500	zmoins = 0.0;
1501	epsiln = sn + ecn;
1502	pc = std::sqrt(std::pow((1.0 + (ecn)/9.3956e2),2) - 1.0) * 9.3956e2;
1503	if(itest == 1) {
1504	G4cout <<"PK::: pc " << pc << G4endl;
1505	}
1506	malpha = 0.0;
1507
1508	// volant:
1509	volant->iv = volant->iv + 1;
1510	volant->acv[volant->iv] = 1.;
1511	volant->zpcv[volant->iv] = 0.;
1512	volant->pcv[volant->iv] = pc;
1513
1514	if(volant->getTotalMass() > 209 && verboseLevel > 0) {
1515	volant->dump();
1516	G4cout <<"DEBUGA Total = " << volant->getTotalMass() << G4endl;
1517	}
1518	}
1519	else {
1520	// fission
1521	// in case of fission-events the fragment nucleus is the mother nucleus
1522	// before fission occurs with excitation energy above the fis.- barrier.
1523	// fission fragment mass distribution is calulated in subroutine fisdis
1524	if (itest == 1) {
1525	G4cout <<"PK::: < fission >" << G4endl;
1526	}
1527	amoins = 0.0;
1528	zmoins = 0.0;
1529	epsiln = ef;
1530
1531	malpha = 0.0;
1532	pc = 0.0;
1533	ff = 1;
1534	// ff = 0; // For testing, allows to disable fission!
1535	}
1536
1537	if (itest == 1) {
1538	G4cout << std::setprecision(9) <<"PK::: SN,SBP,SBA,EF " << sn << " " << sbp << " " << sba <<" " << ef << G4endl;
1539	G4cout << std::setprecision(9) <<"PK::: PROBN,PROBP,PROBA,PROBF,PTOTL " <<" "<< probn <<" "<< probp <<" "<< proba <<" "<< probf <<" "<< ptotl << G4endl;
1540	}
1541
1542	// calculation of the daughter nucleus
1543	af = af - amoins;
1544	zf = zf - zmoins;
1545	ee = ee - epsiln;
1546	if (ee <= 0.01) {
1547	ee = 0.01;
1548	}
1549	mtota = mtota + malpha;
1550
1551	if(ff == 0) {
1552	rnd = randomGenerator->getRandom();
1553	// standardRandom(&rnd,&(hazard->igraine[8]));
1554	ctet1 = 2.0*rnd - 1.0;
1555	rnd = randomGenerator->getRandom();
1556	// standardRandom(&rnd,&(hazard->igraine[4]));
1557	phi1 = rnd2.03.141592654;
1558	stet1 = std::sqrt(1.0 - std::pow(ctet1,2));
1559	volant->xcv[volant->iv] = stet1*std::cos(phi1);
1560	volant->ycv[volant->iv] = stet1*std::sin(phi1);
1561	volant->zcv[volant->iv] = ctet1;
1562	pxeva = pxeva - pc * volant->xcv[volant->iv];
1563	pyeva = pyeva - pc * volant->ycv[volant->iv];
1564	pleva = pleva - pc * ctet1;
1565	}
1566
1567	// condition for end of evaporation
1568	if ((af < 2.5) \|\| (ff == 1)) {
1569	goto evapora100;
1570	}
1571	goto evapora10;
1572
1573	evapora100:
1574	(*zf_par) = zf;
1575	(*af_par) = af;
1576	(*mtota_par) = mtota;
1577	(*pleva_par) = pleva;
1578	(*pxeva_par) = pxeva;
1579	(*pyeva_par) = pyeva;
1580	(*ff_par) = ff;
1581	(*inttype_par) = inttype;
1582	(*inum_par) = inum;
1583
1584	return;
1585	}
1586
1587	void G4Abla::direct(G4double zprf, G4double a, G4double ee, G4double jprf,
1588	G4double probp_par, G4double probn_par, G4double *proba_par,
1589	G4double probf_par, G4double ptotl_par, G4double *sn_par,
1590	G4double sbp_par, G4double sba_par, G4double *ecn_par,
1591	G4double ecp_par,G4double eca_par, G4double *bp_par,
1592	G4double *ba_par, G4int inttype, G4int inum, G4int itest)
1593	{
1594	G4int dummy0 = 0;
1595
1596	G4double probp = (*probp_par);
1597	G4double probn = (*probn_par);
1598	G4double proba = (*proba_par);
1599	G4double probf = (*probf_par);
1600	G4double ptotl = (*ptotl_par);
1601	G4double sn = (*sn_par);
1602	G4double sbp = (*sbp_par);
1603	G4double sba = (*sba_par);
1604	G4double ecn = (*ecn_par);
1605	G4double ecp = (*ecp_par);
1606	G4double eca = (*eca_par);
1607	G4double bp = (*bp_par);
1608	G4double ba = (*ba_par);
1609
1610	// CALCULATION OF PARTICLE-EMISSION PROBABILITIES & FISSION /
1611	// BASED ON THE SIMPLIFIED FORMULAS FOR THE DECAY WIDTH BY /
1612	// MORETTO, ROCHESTER MEETING TO AVOID COMPUTING TIME /
1613	// INTENSIVE INTEGRATION OF THE LEVEL DENSITIES /
1614	// USES EFFECTIVE COULOMB BARRIERS AND AN AVERAGE KINETIC ENERGY/
1615	// OF THE EVAPORATED PARTICLES /
1616	// COLLECTIVE ENHANCMENT OF THE LEVEL DENSITY IS INCLUDED /
1617	// DYNAMICAL HINDRANCE OF FISSION IS INCLUDED BY A STEP FUNCTION/
1618	// APPROXIMATION. SEE A.R. JUNGHANS DIPLOMA THESIS /
1619	// SHELL AND PAIRING STRUCTURES IN THE LEVEL DENSITY IS INCLUDED/
1620
1621	// INPUT:
1622	// ZPRF,A,EE CHARGE, MASS, EXCITATION ENERGY OF COMPOUND
1623	// NUCLEUS
1624	// JPRF ROOT-MEAN-SQUARED ANGULAR MOMENTUM
1625
1626	// DEFORMATIONS AND G.S. SHELL EFFECTS
1627	// COMMON /ECLD/ ECGNZ,ECFNZ,VGSLD,ALPHA
1628
1629	// ECGNZ - GROUND STATE SHELL CORR. FRLDM FOR A SPHERICAL G.S.
1630	// ECFNZ - SHELL CORRECTION FOR THE SADDLE POINT (NOW: == 0)
1631	// VGSLD - DIFFERENCE BETWEEN DEFORMED G.S. AND LDM VALUE
1632	// ALPHA - ALPHA GROUND STATE DEFORMATION (THIS IS NOT BETA2!)
1633	// BETA2 = SQRT((4PI)/5) * ALPHA
1634
1635	//OPTIONS AND PARAMETERS FOR FISSION CHANNEL
1636	//COMMON /FISS/ AKAP,BET,HOMEGA,KOEFF,IFIS,
1637	// OPTSHP,OPTXFIS,OPTLES,OPTCOL
1638	//
1639	// AKAP - HBAR*2/(2 MN * R_0**2) = 10 MEV, R_0 = 1.4 FM
1640	// BET - REDUCED NUCLEAR FRICTION COEFFICIENT IN (1021 S-1)
1641	// HOMEGA - CURVATURE OF THE FISSION BARRIER = 1 MEV
1642	// KOEFF - COEFFICIENT FOR THE LD FISSION BARRIER == 1.0
1643	// IFIS - 0/1 FISSION CHANNEL OFF/ON
1644	// OPTSHP - INTEGER SWITCH FOR SHELL CORRECTION IN MASSES/ENERGY
1645	// = 0 NO MICROSCOPIC CORRECTIONS IN MASSES AND ENERGY
1646	// = 1 SHELL , NO PAIRING
1647	// = 2 PAIRING, NO SHELL
1648	// = 3 SHELL AND PAIRING
1649	// OPTCOL - 0/1 COLLECTIVE ENHANCEMENT SWITCHED ON/OFF
1650	// OPTXFIS- 0,1,2 FOR MYERS & SWIATECKI, DAHLINGER, ANDREYEV
1651	// FISSILITY PARAMETER.
1652	// OPTLES - CONSTANT TEMPERATURE LEVEL DENSITY FOR A,Z > TH-224
1653	// OPTCOL - 0/1 COLLECTIVE ENHANCEMENT OFF/ON
1654
1655	// LEVEL DENSITY PARAMETERS
1656	// COMMON /ALD/ AV,AS,AK,OPTAFAN
1657	// AV,AS,AK - VOLUME,SURFACE,CURVATURE DEPENDENCE OF THE
1658	// LEVEL DENSITY PARAMETER
1659	// OPTAFAN - 0/1 AF/AN >=1 OR AF/AN ==1
1660	// RECOMMENDED IS OPTAFAN = 0
1661
1662	// FISSION BARRIERS
1663	// COMMON /FB/ EFA
1664	// EFA - ARRAY OF FISSION BARRIERS
1665
1666
1667	// OUTPUT: PROBN,PROBP,PROBA,PROBF,PTOTL:
1668	// - EMISSION PROBABILITIES FOR N EUTRON, P ROTON, A LPHA
1669	// PARTICLES, F ISSION AND NORMALISATION
1670	// SN,SBP,SBA: SEPARATION ENERGIES N P A
1671	// INCLUDING EFFECTIVE BARRIERS
1672	// ECN,ECP,ECA,BP,BA
1673	// - AVERAGE KINETIC ENERGIES (2*T) AND EFFECTIVE BARRIERS
1674
1675	static G4double bk = 0.0;
1676	static G4int afp = 0;
1677	static G4double at = 0.0;
1678	static G4double bs = 0.0;
1679	static G4double bshell = 0.0;
1680	static G4double cf = 0.0;
1681	static G4double dconst = 0.0;
1682	static G4double defbet = 0.0;
1683	static G4double denomi = 0.0;
1684	static G4double densa = 0.0;
1685	static G4double densf = 0.0;
1686	static G4double densg = 0.0;
1687	static G4double densn = 0.0;
1688	static G4double densp = 0.0;
1689	static G4double edyn = 0.0;
1690	static G4double eer = 0.0;
1691	static G4double ef = 0.0;
1692	static G4double ft = 0.0;
1693	static G4double ga = 0.0;
1694	static G4double gf = 0.0;
1695	static G4double gn = 0.0;
1696	static G4double gngf = 0.0;
1697	static G4double gp = 0.0;
1698	static G4double gsum = 0.0;
1699	static G4double hbar = 6.582122e-22; // = 0.0;
1700	static G4double iflag = 0.0;
1701	static G4int il = 0;
1702	static G4int imaxwell = 0;
1703	static G4int in = 0;
1704	static G4int iz = 0;
1705	static G4int j = 0;
1706	static G4int k = 0;
1707	static G4double ma1z = 0.0;
1708	static G4double ma1z1 = 0.0;
1709	static G4double ma4z2 = 0.0;
1710	static G4double maz = 0.0;
1711	static G4double nprf = 0.0;
1712	static G4double nt = 0.0;
1713	static G4double parc = 0.0;
1714	static G4double pi = 3.14159265;
1715	static G4double pt = 0.0;
1716	static G4double ra = 0.0;
1717	static G4double rat = 0.0;
1718	static G4double refmod = 0.0;
1719	static G4double rf = 0.0;
1720	static G4double rn = 0.0;
1721	static G4double rnd = 0.0;
1722	static G4double rnt = 0.0;
1723	static G4double rp = 0.0;
1724	static G4double rpt = 0.0;
1725	static G4double sa = 0.0;
1726	static G4double sbf = 0.0;
1727	static G4double sbfis = 0.0;
1728	static G4double segs = 0.0;
1729	static G4double selmax = 0.0;
1730	static G4double sp = 0.0;
1731	static G4double tauc = 0.0;
1732	static G4double tconst = 0.0;
1733	static G4double temp = 0.0;
1734	static G4double ts1 = 0.0;
1735	static G4double tsum = 0.0;
1736	static G4double wf = 0.0;
1737	static G4double wfex = 0.0;
1738	static G4double xx = 0.0;
1739	static G4double y = 0.0;
1740
1741	imaxwell = 1;
1742	inttype = 0;
1743
1744	// limiting of excitation energy where fission occurs
1745	// Note, this is not the dynamical hindrance (see end of routine)
1746	edyn = 1000.0;
1747
1748	// no limit if statistical model is calculated.
1749	if (fiss->bet <= 1.0e-16) {
1750	edyn = 10000.0;
1751	}
1752
1753	// just a change of name until the end of this subroutine
1754	eer = ee;
1755	if(verboseLevel > 2)
1756	G4cout << __FILE__ << ":" << __LINE__ << " eer = " << eer << G4endl;
1757	if (inum == 1) {
1758	ilast = 1;
1759	}
1760
1761	// calculation of masses
1762	// refmod = 1 ==> myers,swiatecki model
1763	// refmod = 0 ==> weizsaecker model
1764	refmod = 1; // Default = 1
1765
1766	if (refmod == 1) {
1767	mglms(a,zprf,fiss->optshp,&maz);
1768	mglms(a-1.0,zprf,fiss->optshp,&ma1z);
1769	mglms(a-1.0,zprf-1.0,fiss->optshp,&ma1z1);
1770	mglms(a-4.0,zprf-2.0,fiss->optshp,&ma4z2);
1771	}
1772	else {
1773	mglw(a,zprf,&maz);
1774	mglw(a-1.0,zprf,&ma1z);
1775	mglw(a-1.0,zprf-1.0,&ma1z1);
1776	mglw(a-4.0,zprf-2.0,&ma4z2);
1777	}
1778
1779	// separation energies and effective barriers
1780	sn = ma1z - maz;
1781	sp = ma1z1 - maz;
1782	sa = ma4z2 - maz - 28.29688;
1783	if (zprf < 1.0e0) {
1784	sbp = 1.0e75;
1785	goto direct30;
1786	}
1787
1788	// parameterisation gaimard:
1789	// bp = 1.44(zprf-1.d0)/(1.22std::pow((a - 1.0),(1.0/3.0))+5.6)
1790	// parameterisation khs (12-99)
1791	bp = 1.44(zprf - 1.0)/(2.1std::pow((a - 1.0),(1.0/3.0)) + 0.0);
1792
1793	sbp = sp + bp;
1794	if (a-4.0 <= 0.0) {
1795	sba = 1.0e+75;
1796	goto direct30;
1797	}
1798
1799	// new effective barrier for alpha evaporation d=6.1: khs
1800	// ba = 2.88d0(zprf-2.d0)/(1.22d0(a-4.d0)**(1.d0/3.d0)+6.1d0)
1801	// parametrisation khs (12-99)
1802	ba = 2.88(zprf - 2.0)/(2.2std::pow((a - 4.0),(1.0/3.0)) + 0.0);
1803
1804	sba = sa + ba;
1805	direct30:
1806
1807	// calculation of surface and curvature integrals needed to
1808	// to calculate the level density parameter (in densniv)
1809	if (fiss->ifis > 0) {
1810	k = idnint(zprf);
1811	j = idnint(a - zprf);
1812
1813	// now ef is calculated from efa that depends on the subroutine
1814	// barfit which takes into account the modification on the ang. mom.
1815	// jb mvr 6-aug-1999
1816	// note *** shell correction! (ecgnz) jb mvr 20-7-1999
1817	il = idnint(jprf);
1818	barfit(k,k+j,il,&sbfis,&segs,&selmax);
1819	if ((fiss->optshp == 1) \|\| (fiss->optshp == 3)) {
1820	// fb->efa[k][j] = G4double(sbfis) - ecld->ecgnz[j][k];
1821	// fb->efa[j][k] = G4double(sbfis) - ecld->ecgnz[j][k];
1822	fb->efa[j][k] = double(sbfis) - ecld->ecgnz[j][k];
1823	}
1824	else {
1825	// fb->efa[k][j] = G4double(sbfis);
1826	fb->efa[j][k] = double(sbfis);
1827	}
1828	// ef = fb->efa[k][j];
1829	ef = fb->efa[j][k];
1830
1831	// to avoid negative values for impossible nuclei
1832	// the fission barrier is set to zero if smaller than zero.
1833	// if (fb->efa[k][j] < 0.0) {
1834	// fb->efa[k][j] = 0.0;
1835	// }
1836	if (fb->efa[j][k] < 0.0) {
1837	if(verboseLevel > 2) {
1838	G4cout <<"Setting fission barrier to 0" << G4endl;
1839	}
1840	fb->efa[j][k] = 0.0;
1841	}
1842
1843	// factor with jprf should be 0.0025d0 - 0.01d0 for
1844	// approximate influence of ang. momentum on bfis a.j. 22.07.96
1845	// 0.0 means no angular momentum
1846
1847	if (ef < 0.0) {
1848	ef = 0.0;
1849	}
1850	xx = fissility((k+j),k,fiss->optxfis);
1851
1852	y = 1.00 - xx;
1853	if (y < 0.0) {
1854	y = 0.0;
1855	}
1856	if (y > 1.0) {
1857	y = 1.0;
1858	}
1859	bs = bipol(1,y);
1860	bk = bipol(2,y);
1861	}
1862	else {
1863	ef = 1.0e40;
1864	bs = 1.0;
1865	bk = 1.0;
1866	}
1867	sbf = ee - ef;
1868
1869	afp = idnint(a);
1870	iz = idnint(zprf);
1871	in = afp - iz;
1872	bshell = ecld->ecfnz[in][iz];
1873
1874	// ld saddle point deformation
1875	// here: beta2 = std::sqrt(5/(4pi)) * alpha2
1876
1877	// for the ground state def. 1.5d0 should be used
1878	// because this was just the factor to produce the
1879	// alpha-deformation table 1.5d0 should be used
1880	// a.r.j. 6.8.97
1881	defbet = 1.58533e0 * spdef(idnint(a),idnint(zprf),fiss->optxfis);
1882
1883	// level density and temperature at the saddle point
1884	// G4cout <<"a = " << a << G4endl;
1885	// G4cout <<"zprf = " << zprf << G4endl;
1886	// G4cout <<"ee = " << ee << G4endl;
1887	// G4cout <<"ef = " << ef << G4endl;
1888	// G4cout <<"bshell = " << bshell << G4endl;
1889	// G4cout <<"bs = " << bs << G4endl;
1890	// G4cout <<"bk = " << bk << G4endl;
1891	// G4cout <<"defbet = " << defbet << G4endl;
1892	densniv(a,zprf,ee,ef,&densf,bshell,bs,bk,&temp,int(fiss->optshp),int(fiss->optcol),defbet);
1893	// G4cout <<"densf = " << densf << G4endl;
1894	// G4cout <<"temp = " << temp << G4endl;
1895	ft = temp;
1896	if (iz >= 2) {
1897	bshell = ecld->ecgnz[in][iz-1] - ecld->vgsld[in][iz-1];
1898	defbet = 1.5 * (ecld->alpha[in][iz-1]);
1899
1900	// level density and temperature in the proton daughter
1901	densniv(a-1.0,zprf-1.0e0,ee,sbp,&densp, bshell,1.e0,1.e0,&temp,int(fiss->optshp),int(fiss->optcol),defbet);
1902	pt = temp;
1903	if (imaxwell == 1) {
1904	// valentina - random kinetic energy in a maxwelliam distribution
1905	// modif juin/2002 a.b. c.v. for light targets; limit on the energy
1906	// from the maxwell distribution.
1907	rpt = pt;
1908	ecp = 2.0 * pt;
1909	if(rpt <= 1.0e-3) {
1910	goto direct2914;
1911	}
1912	iflag = 0;
1913	direct1914:
1914	ecp = fmaxhaz(rpt);
1915	iflag = iflag + 1;
1916	if(iflag >= 10) {
1917	rnd = randomGenerator->getRandom();
1918	// standardRandom(&rnd,&(hazard->igraine[5]));
1919	ecp=std::sqrt(rnd)*(eer-sbp);
1920	goto direct2914;
1921	}
1922	if((ecp+sbp) > eer) {
1923	goto direct1914;
1924	}
1925	}
1926	else {
1927	ecp = 2.0 * pt;
1928	}
1929
1930	direct2914:
1931	dummy0 = 0;
1932	// G4cout <<""<<G4endl;
1933	}
1934	else {
1935	densp = 0.0;
1936	ecp = 0.0;
1937	pt = 0.0;
1938	}
1939
1940	if (in >= 2) {
1941	bshell = ecld->ecgnz[in-1][iz] - ecld->vgsld[in-1][iz];
1942	defbet = 1.5e0 * (ecld->alpha[in-1][iz]);
1943
1944	// level density and temperature in the neutron daughter
1945	densniv(a-1.0,zprf,ee,sn,&densn,bshell, 1.e0,1.e0,&temp,int(fiss->optshp),int(fiss->optcol),defbet);
1946	nt = temp;
1947	if(itest == 1)
1948	G4cout <<"PK::: nt = " << nt <<G4endl;
1949	if (imaxwell == 1) {
1950	// valentina - random kinetic energy in a maxwelliam distribution
1951	// modif juin/2002 a.b. c.v. for light targets; limit on the energy
1952	// from the maxwell distribution.
1953	rnt = nt;
1954	ecn = 2.0 * nt;
1955	if(rnt <= 1.e-3) {
1956	goto direct2915;
1957	}
1958
1959	iflag=0;
1960	direct1915:
1961	ecn = fmaxhaz(rnt);
1962	if(verboseLevel > 2) {
1963	G4cout <<"rnt = " << rnt << G4endl;
1964	G4cout << __FILE__ << ":" << __LINE__ << " ecn = " << ecn << G4endl;
1965	}
1966	iflag=iflag+1;
1967	if(iflag >= 10) {
1968	rnd = randomGenerator->getRandom();
1969	// standardRandom(&rnd,&(hazard->igraine[6]));
1970	ecn = std::sqrt(rnd)*(eer-sn);
1971	if(verboseLevel > 2)
1972	G4cout << __FILE__ << ":" << __LINE__ << " ecn = " << ecn << G4endl;
1973	goto direct2915;
1974	}
1975	if((ecn+sn) > eer) {
1976	goto direct1915;
1977	}
1978	}
1979	else {
1980	ecn = 2.e0 * nt;
1981	if(verboseLevel > 2)
1982	G4cout << __FILE__ << ":" << __LINE__ << " ecn = " << ecn << G4endl;
1983	}
1984	// if((ecn + sn) <= eer) {
1985	// ecn = 2.0 * nt;
1986	// G4cout << __FILE__ << ":" << __LINE__ << " ecn = " << ecn << G4endl;
1987	// }
1988	direct2915:
1989	dummy0 = 0;
1990	// G4cout <<"" <<G4endl;
1991	}
1992	else {
1993	densn = 0.0;
1994	ecn = 0.0;
1995	nt = 0.0;
1996	}
1997
1998	if ((in >= 3) && (iz >= 3)) {
1999	bshell = ecld->ecgnz[in-2][iz-2] - ecld->vgsld[in-2][iz-2];
2000	defbet = 1.5 * (ecld->alpha[in-2][iz-2]);
2001
2002	// For debugging:
2003	// bshell = -10.7;
2004	// defbet = -0.06105;
2005	// G4cout <<"ecgnz N = " << in-2 << G4endl;
2006	// G4cout <<"ecgnz Z = " << iz-2 << G4endl;
2007	// G4cout <<"bshell = " << bshell << G4endl;
2008	// G4cout <<"defbet = " << defbet << G4endl;
2009	// level density and temperature in the alpha daughter
2010	densniv(a-4.0,zprf-2.0e0,ee,sba,&densa,bshell,1.e0,1.e0,&temp,int(fiss->optshp),int(fiss->optcol),defbet);
2011	// G4cout <<"densa = " << densa << G4endl;
2012	// G4cout <<"temp = " << temp << G4endl;
2013
2014	// valentina - random kinetic energy in a maxwelliam distribution
2015	at = temp;
2016	if (imaxwell == 1) {
2017	// modif juin/2002 a.b. c.v. for light targets; limit on the energy
2018	// from the maxwell distribution.
2019	rat = at;
2020	eca= 2.e0 * at;
2021	if(rat <= 1.e-3) {
2022	goto direct2916;
2023	}
2024	iflag=0;
2025	direct1916:
2026	eca = fmaxhaz(rat);
2027	iflag=iflag+1;
2028	if(iflag >= 10) {
2029	rnd = randomGenerator->getRandom();
2030	// standardRandom(&rnd,&(hazard->igraine[7]));
2031	eca=std::sqrt(rnd)*(eer-sba);
2032	goto direct2916;
2033	}
2034	if((eca+sba) > eer) {
2035	goto direct1916;
2036	}
2037	}
2038	else {
2039	eca = 2.0 * at;
2040	}
2041	direct2916:
2042	dummy0 = 0;
2043	// G4cout <<"" << G4endl;
2044	}
2045	else {
2046	densa = 0.0;
2047	eca = 0.0;
2048	at = 0.0;
2049	}
2050	//} // PK
2051
2052	// special treatment for unbound nuclei
2053	if (sn < 0.0) {
2054	probn = 1.0;
2055	probp = 0.0;
2056	proba = 0.0;
2057	probf = 0.0;
2058	goto direct70;
2059	}
2060	if (sbp < 0.0) {
2061	probp = 1.0;
2062	probn = 0.0;
2063	proba = 0.0;
2064	probf = 0.0;
2065	goto direct70;
2066	}
2067
2068	if ((a < 50.e0) \|\| (ee > edyn)) { // no fission if e*> edyn or mass < 50
2069	// G4cout <<"densf = 0.0" << G4endl;
2070	densf = 0.e0;
2071	}
2072
2073	bshell = ecld->ecgnz[in][iz] - ecld->vgsld[in][iz];
2074	defbet = 1.5e0 * (ecld->alpha[in][iz]);
2075
2076	// compound nucleus level density
2077	densniv(a,zprf,ee,0.0e0,&densg,bshell,1.e0,1.e0,&temp,int(fiss->optshp),int(fiss->optcol),defbet);
2078
2079	if ( densg > 0.e0) {
2080	// calculation of the partial decay width
2081	// used for both the time scale and the evaporation decay width
2082	gp = (std::pow(a,(2.0/3.0))/fiss->akap)densp/densg/pistd::pow(pt,2);
2083	gn = (std::pow(a,(2.0/3.0))/fiss->akap)densn/densg/pistd::pow(nt,2);
2084	ga = (std::pow(a,(2.0/3.0))/fiss->akap)densa/densg/pi2.0*std::pow(at,2);
2085	gf = densf/densg/pi/2.0*ft;
2086
2087	if(itest == 1) {
2088	G4cout <<"gn,gp,ga,gf " << gn <<","<< gp <<","<< ga <<","<< gf << G4endl;
2089	}
2090	}
2091	else {
2092	if(verboseLevel > 2) {
2093	G4cout <<"direct: densg <= 0.e0 " << a <<","<< zprf <<","<< ee << G4endl;
2094	}
2095	}
2096
2097	gsum = ga + gp + gn;
2098	if (gsum > 0.0) {
2099	ts1 = hbar / gsum;
2100	}
2101	else {
2102	ts1 = 1.0e99;
2103	}
2104
2105	if (inum > ilast) { // new event means reset the time scale
2106	tsum = 0;
2107	}
2108
2109	// calculate the relative probabilities for all decay channels
2110	if (densf == 0.0) {
2111	if (densp == 0.0) {
2112	if (densn == 0.0) {
2113	if (densa == 0.0) {
2114	// no reaction is possible
2115	probf = 0.0;
2116	probp = 0.0;
2117	probn = 0.0;
2118	proba = 0.0;
2119	goto direct70;
2120	}
2121
2122	// alpha evaporation is the only open channel
2123	rf = 0.0;
2124	rp = 0.0;
2125	rn = 0.0;
2126	ra = 1.0;
2127	goto direct50;
2128	}
2129
2130	// alpha emission and neutron emission
2131	rf = 0.0;
2132	rp = 0.0;
2133	rn = 1.0;
2134	ra = densa2.0/densnstd::pow((at/nt),2);
2135	goto direct50;
2136	}
2137	// alpha, proton and neutron emission
2138	rf = 0.0;
2139	rp = 1.0;
2140	rn = densn/densp*std::pow((nt/pt),2);
2141	ra = densa2.0/denspstd::pow((at/pt),2);
2142	goto direct50;
2143	}
2144
2145	// here fission has taken place
2146	rf = 1.0;
2147
2148	// cramers and weidenmueller factors for the dynamical hindrances of
2149	// fission
2150	if (fiss->bet <= 1.0e-16) {
2151	cf = 1.0;
2152	wf = 1.0;
2153	}
2154	else if (sbf > 0.0e0) {
2155	cf = cram(fiss->bet,fiss->homega);
2156	// if fission barrier ef=0.d0 then fission is the only possible
2157	// channel. to avoid std::log(0) in function tau
2158	// a.j. 7/28/93
2159	if (ef <= 0.0) {
2160	rp = 0.0;
2161	rn = 0.0;
2162	ra = 0.0;
2163	goto direct50;
2164	}
2165	else {
2166	// transient time tau()
2167	tauc = tau(fiss->bet,fiss->homega,ef,ft);
2168	}
2169	wfex = (tauc - tsum)/ts1;
2170
2171	if (wfex < 0.0) {
2172	wf = 1.0;
2173	}
2174	else {
2175	wf = std::exp( -wfex);
2176	}
2177	}
2178	else {
2179	cf=1.0;
2180	wf=1.0;
2181	}
2182
2183	if(verboseLevel > 2) {
2184	G4cout <<"tsum,wf,cf " << tsum <<","<< wf <<","<< cf << G4endl;
2185	}
2186
2187	tsum = tsum + ts1;
2188
2189	// change by g.k. and a.h. 5.9.95
2190	tconst = 0.7;
2191	dconst = 12.0/std::sqrt(a);
2192	nprf = a - zprf;
2193
2194	if (fiss->optshp >= 2) { //then
2195	parite(nprf,&parc);
2196	dconst = dconst*parc;
2197	}
2198	else {
2199	dconst= 0.0;
2200	}
2201	if ((ee <= 17.e0) && (fiss->optles == 1) && (iz >= 90) && (in >= 134)) { //then
2202	// constant changed to 5.0 accord to moretto & vandenbosch a.j. 19.3.96
2203	gngf = std::pow(a,(2.0/3.0))tconst/10.0std::exp((ef-sn+dconst)/tconst);
2204
2205	// if the excitation energy is so low that densn=0 ==> gn = 0
2206	// fission remains the only channel.
2207	// a. j. 10.1.94
2208	if (gn == 0.0) {
2209	rn = 0.0;
2210	rp = 0.0;
2211	ra = 0.0;
2212	}
2213	else {
2214	rn=gngf;
2215	rp=gngf*gp/gn;
2216	ra=gngf*ga/gn;
2217	}
2218	} else {
2219	rn = gn/(gf*cf);
2220	// G4cout <<"rn = " << G4endl;
2221	// G4cout <<"gn = " << gn << " gf = " << gf << " cf = " << cf << G4endl;
2222	rp = gp/(gf*cf);
2223	ra = ga/(gf*cf);
2224	}
2225	direct50:
2226	// relative decay probabilities
2227	denomi = rp+rn+ra+rf;
2228	// decay probabilities after transient time
2229	probf = rf/denomi;
2230	probp = rp/denomi;
2231	probn = rn/denomi;
2232	proba = ra/denomi;
2233
2234	// new treatment of grange-weidenmueller factor, 5.1.2000, khs !!!
2235
2236	// decay probabilites with transient time included
2237	probf = probf * wf;
2238	if(probf == 1.0) {
2239	probp = 0.0;
2240	probn = 0.0;
2241	proba = 0.0;
2242	}
2243	else {
2244	probp = probp * (wf + (1.e0-wf)/(1.e0-probf));
2245	probn = probn * (wf + (1.e0-wf)/(1.e0-probf));
2246	proba = proba * (wf + (1.e0-wf)/(1.e0-probf));
2247	}
2248	direct70:
2249	ptotl = probp+probn+proba+probf;
2250
2251	ee = eer;
2252	ilast = inum;
2253
2254	// Return values:
2255	(*probp_par) = probp;
2256	(*probn_par) = probn;
2257	(*proba_par) = proba;
2258	(*probf_par) = probf;
2259	(*ptotl_par) = ptotl;
2260	(*sn_par) = sn;
2261	(*sbp_par) = sbp;
2262	(*sba_par) = sba;
2263	(*ecn_par) = ecn;
2264	(*ecp_par) = ecp;
2265	(*eca_par) = eca;
2266	(*bp_par) = bp;
2267	(*ba_par) = ba;
2268	}
2269
2270	void G4Abla::densniv(G4double a, G4double z, G4double ee, G4double esous, G4double *dens, G4double bshell, G4double bs, G4double bk,
2271	G4double *temp, G4int optshp, G4int optcol, G4double defbet)
2272	{
2273	// 1498 C
2274	// 1499 C INPUT:
2275	// 1500 C A,EE,ESOUS,OPTSHP,BS,BK,BSHELL,DEFBET
2276	// 1501 C
2277	// 1502 C LEVEL DENSITY PARAMETERS
2278	// 1503 C COMMON /ALD/ AV,AS,AK,OPTAFAN
2279	// 1504 C AV,AS,AK - VOLUME,SURFACE,CURVATURE DEPENDENCE OF THE
2280	// 1505 C LEVEL DENSITY PARAMETER
2281	// 1506 C OPTAFAN - 0/1 AF/AN >=1 OR AF/AN ==1
2282	// 1507 C RECOMMENDED IS OPTAFAN = 0
2283	// 1508 C---------------------------------------------------------------------
2284	// 1509 C OUTPUT: DENS,TEMP
2285	// 1510 C
2286	// 1511 C ____________________________________________________________________
2287	// 1512 C /
2288	// 1513 C / PROCEDURE FOR CALCULATING THE STATE DENSITY OF A COMPOUND NUCLEUS
2289	// 1514 C /____________________________________________________________________
2290	// 1515 C
2291	// 1516 INTEGER AFP,IZ,OPTSHP,OPTCOL,J,OPTAFAN
2292	// 1517 REAL*8 A,EE,ESOUS,DENS,E,Y0,Y1,Y2,Y01,Y11,Y21,PA,BS,BK,TEMP
2293	// 1518 C=====INSERTED BY KUDYAEV===============================================
2294	// 1519 COMMON /ALD/ AV,AS,AK,OPTAFAN
2295	// 1520 REAL*8 ECR,ER,DELTAU,Z,DELTPP,PARA,PARZ,FE,HE,ECOR,ECOR1,Pi6
2296	// 1521 REAL*8 BSHELL,DELTA0,AV,AK,AS,PONNIV,PONFE,DEFBET,QR,SIG,FP
2297	// 1522 C=======================================================================
2298	// 1523 C
2299	// 1524 C
2300	// 1525 C-----------------------------------------------------------------------
2301	// 1526 C A MASS NUMBER OF THE DAUGHTER NUCLEUS
2302	// 1527 C EE EXCITATION ENERGY OF THE MOTHER NUCLEUS
2303	// 1528 C ESOUS SEPARATION ENERGY PLUS EFFECTIVE COULOMB BARRIER
2304	// 1529 C DENS STATE DENSITY OF DAUGHTER NUCLEUS AT EE-ESOUS-EC
2305	// 1530 C BSHELL SHELL CORRECTION
2306	// 1531 C TEMP NUCLEAR TEMPERATURE
2307	// 1532 C E LOCAL EXCITATION ENERGY OF THE DAUGHTER NUCLEUS
2308	// 1533 C E1 LOCAL HELP VARIABLE
2309	// 1534 C Y0,Y1,Y2,Y01,Y11,Y21
2310	// 1535 C LOCAL HELP VARIABLES
2311	// 1536 C PA LOCAL STATE-DENSITY PARAMETER
2312	// 1537 C EC KINETIC ENERGY OF EMITTED PARTICLE WITHOUT
2313	// 1538 C COULOMB REPULSION
2314	// 1539 C IDEN FAKTOR FOR SUBSTRACTING KINETIC ENERGY IDEN*TEMP
2315	// 1540 C DELTA0 PAIRING GAP 12 FOR GROUND STATE
2316	// 1541 C 14 FOR SADDLE POINT
2317	// 1542 C EITERA HELP VARIABLE FOR TEMPERATURE ITERATION
2318	// 1543 C-----------------------------------------------------------------------
2319	// 1544 C
2320	// 1545 C
2321	G4double afp = 0.0;
2322	G4double delta0 = 0.0;
2323	G4double deltau = 0.0;
2324	G4double deltpp = 0.0;
2325	G4double e = 0.0;
2326	G4double ecor = 0.0;
2327	G4double ecor1 = 0.0;
2328	G4double ecr = 0.0;
2329	G4double er = 0.0;
2330	G4double fe = 0.0;
2331	G4double fp = 0.0;
2332	G4double he = 0.0;
2333	G4double iz = 0.0;
2334	G4double pa = 0.0;
2335	G4double para = 0.0;
2336	G4double parz = 0.0;
2337	G4double ponfe = 0.0;
2338	G4double ponniv = 0.0;
2339	G4double qr = 0.0;
2340	G4double sig = 0.0;
2341	G4double y01 = 0.0;
2342	G4double y11 = 0.0;
2343	G4double y2 = 0.0;
2344	G4double y21 = 0.0;
2345	G4double y1 = 0.0;
2346	G4double y0 = 0.0;
2347
2348	G4double pi6 = std::pow(3.1415926535,2) / 6.0;
2349	ecr=10.0;
2350	er=28.0;
2351	afp=idnint(a);
2352	iz=idnint(z);
2353
2354	// level density parameter
2355	if((ald->optafan == 1)) {
2356	pa = (ald->av)a + (ald->as)std::pow(a,(2.e0/3.e0)) + (ald->ak)*std::pow(a,(1.e0/3.e0));
2357	}
2358	else {
2359	pa = (ald->av)a + (ald->as)bsstd::pow(a,(2.e0/3.e0)) + (ald->ak)bk*std::pow(a,(1.e0/3.e0));
2360	}
2361
2362	fp = 0.01377937231e0 * std::pow(a,(5.e0/3.e0)) * (1.e0 + defbet/3.e0);
2363
2364	// pairing corrections
2365	if (bs > 1.0) {
2366	delta0 = 14;
2367	}
2368	else {
2369	delta0 = 12;
2370	}
2371
2372	if (esous > 1.0e30) {
2373	(*dens) = 0.0;
2374	(*temp) = 0.0;
2375	goto densniv100;
2376	}
2377
2378	e = ee - esous;
2379
2380	if (e < 0.0) {
2381	(*dens) = 0.0;
2382	(*temp) = 0.0;
2383	goto densniv100;
2384	}
2385
2386	// shell corrections
2387	if (optshp > 0) {
2388	deltau = bshell;
2389	if (optshp == 2) {
2390	deltau = 0.0;
2391	}
2392	if (optshp >= 2) {
2393	// pairing energy shift with condensation energy a.r.j. 10.03.97
2394	// deltpp = -0.25e0* (delta0/std::pow(std::sqrt(a),2)) * pa /pi6 + 2.e0*delta0/std::sqrt(a);
2395	deltpp = -0.25e0* std::pow((delta0/std::sqrt(a)),2) * pa /pi6 + 2.e0*delta0/std::sqrt(a);
2396
2397	parite(a,&para);
2398	if (para < 0.0) {
2399	e = e - delta0/std::sqrt(a);
2400	} else {
2401	parite(z,&parz);
2402	if (parz > 0.e0) {
2403	e = e - 2.0*delta0/std::sqrt(a);
2404	} else {
2405	e = e;
2406	}
2407	}
2408	} else {
2409	deltpp = 0.0;
2410	}
2411	} else {
2412	deltau = 0.0;
2413	deltpp = 0.0;
2414	}
2415	if (e < 0.0) {
2416	e = 0.0;
2417	(*temp) = 0.0;
2418	}
2419
2420	// washing out is made stronger ! g.k. 3.7.96
2421	ponfe = -2.5pae*std::pow(a,(-4.0/3.0));
2422
2423	if (ponfe < -700.0) {
2424	ponfe = -700.0;
2425	}
2426	fe = 1.0 - std::exp(ponfe);
2427	if (e < ecr) {
2428	// priv. comm. k.-h. schmidt
2429	he = 1.0 - std::pow((1.0 - e/ecr),2);
2430	}
2431	else {
2432	he = 1.0;
2433	}
2434
2435	// Excitation energy corrected for pairing and shell effects
2436	// washing out with excitation energy is included.
2437	ecor = e + deltaufe + deltpphe;
2438
2439	if (ecor <= 0.1) {
2440	ecor = 0.1;
2441	}
2442
2443	// statt 170.d0 a.r.j. 8.11.97
2444
2445	// iterative procedure according to grossjean and feldmeier
2446	// to avoid the singularity e = 0
2447	if (ee < 5.0) {
2448	y1 = std::sqrt(pa*ecor);
2449	for(int j = 0; j < 5; j++) {
2450	y2 = paecor(1.e0-std::exp(-y1));
2451	y1 = std::sqrt(y2);
2452	}
2453
2454	y0 = pa/y1;
2455	(*temp)=1.0/y0;
2456	(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;
2457	if (ecor < 1.0) {
2458	ecor1=1.0;
2459	y11 = std::sqrt(pa*ecor1);
2460	for(int j = 0; j < 7; j++) {
2461	y21 = paecor1(1.0-std::exp(-y11));
2462	y11 = std::sqrt(y21);
2463	}
2464
2465	y01 = pa/y11;
2466	(dens) = (dens)*std::pow((y01/y0),1.5);
2467	(temp) = (temp)*std::pow((y01/y0),1.5);
2468	}
2469	}
2470	else {
2471	ponniv = 2.0std::sqrt(paecor);
2472	if (ponniv > 700.0) {
2473	ponniv = 700.0;
2474	}
2475
2476	// fermi gas state density
2477	(dens) = std::pow(pa,(-0.25e0))std::pow(ecor,(-1.25e0))std::exp(ponniv) 0.1477045e0;
2478	(*temp) = std::sqrt(ecor/pa);
2479	}
2480	densniv100:
2481
2482	// spin cutoff parameter
2483	sig = fp * (*temp);
2484
2485	// collective enhancement
2486	if (optcol == 1) {
2487	qrot(z,a,defbet,sig,ecor,&qr);
2488	}
2489	else {
2490	qr = 1.0;
2491	}
2492
2493	(dens) = (dens) * qr;
2494	if(verboseLevel > 2) {
2495	G4cout <<"PK::: dens = " << (*dens) << G4endl;
2496	G4cout <<"PK::: AFP, IZ, ECOR, ECOR1 " << afp << " " << iz << " " << ecor << " " << ecor1 << G4endl;
2497	}
2498	}
2499
2500
2501	G4double G4Abla::bfms67(G4double zms, G4double ams)
2502	{
2503	// This subroutine calculates the fission barriers
2504	// of the liquid-drop model of Myers and Swiatecki (1967).
2505	// Analytic parameterization of Dahlinger 1982
2506	// replaces tables. Barrier heights from Myers and Swiatecki !!!
2507
2508	G4double nms = 0.0, ims = 0.0, ksims = 0.0, xms = 0.0, ums = 0.0;
2509
2510	nms = ams - zms;
2511	ims = (nms-zms)/ams;
2512	ksims= 50.15e0 * (1.- 1.78 * std::pow(ims,2));
2513	xms = std::pow(zms,2) / (ams * ksims);
2514	ums = 0.368e0-5.057e0xms+8.93e0std::pow(xms,2)-8.71*std::pow(xms,3);
2515	return(0.7322e0std::pow(zms,2)/std::pow(ams,(0.333333e0))std::pow(10.e0,ums));
2516	}
2517
2518	void G4Abla::lpoly(G4double x, G4int n, G4double pl[])
2519	{
2520	// THIS SUBROUTINE CALCULATES THE ORDINARY LEGENDRE POLYNOMIALS OF
2521	// ORDER 0 TO N-1 OF ARGUMENT X AND STORES THEM IN THE VECTOR PL.
2522	// THEY ARE CALCULATED BY RECURSION RELATION FROM THE FIRST TWO
2523	// POLYNOMIALS.
2524	// WRITTEN BY A.J.SIERK LANL T-9 FEBRUARY, 1984
2525
2526	// NOTE: PL AND X MUST BE DOUBLE PRECISION ON 32-BIT COMPUTERS!
2527
2528	pl[0] = 1.0;
2529	pl[1] = x;
2530
2531	for(int i = 2; i < n; i++) {
2532	pl[i] = ((2double(i+1) - 3.0)xpl[i-1] - (double(i+1) - 2.0)pl[i-2])/(double(i+1)-1.0);
2533	}
2534	}
2535
2536	G4double G4Abla::eflmac(G4int ia, G4int iz, G4int flag, G4int optshp)
2537	{
2538	// CHANGED TO CALCULATE TOTAL BINDING ENERGY INSTEAD OF MASS EXCESS.
2539	// SWITCH FOR PAIRING INCLUDED AS WELL.
2540	// BINDING = EFLMAC(IA,IZ,0,OPTSHP)
2541	// FORTRAN TRANSCRIPT OF /U/GREWE/LANG/EEX/FRLDM.C
2542	// A.J. 15.07.96
2543
2544	// this function will calculate the liquid-drop nuclear mass for spheri
2545	// configuration according to the preprint NUCLEAR GROUND-STATE
2546	// MASSES and DEFORMATIONS by P. M"oller et al. from August 16, 1993 p.
2547	// All constants are taken from this publication for consistency.
2548
2549	// Parameters:
2550	// a: nuclear mass number
2551	// z: nuclear charge
2552	// flag: 0 - return mass excess
2553	// otherwise - return pairing (= -1/2 dpn + 1/2 (Dp + Dn))
2554
2555	G4double eflmacResult = 0.0;
2556
2557	G4int in = 0;
2558	G4double z = 0.0, n = 0.0, a = 0.0, av = 0.0, as = 0.0;
2559	G4double a0 = 0.0, c1 = 0.0, c4 = 0.0, b1 = 0.0, b3 = 0.0;
2560	G4double f = 0.0, ca = 0.0, w = 0.0, dp = 0.0, dn = 0.0, dpn = 0.0, efl = 0.0;
2561	G4double rmac = 0.0, bs = 0.0, h = 0.0, r0 = 0.0, kf = 0.0, ks = 0.0;
2562	G4double kv = 0.0, rp = 0.0, ay = 0.0, aden = 0.0, x0 = 0.0, y0 = 0.0;
2563	G4double mh = 0.0, mn = 0.0, esq = 0.0, ael = 0.0, i = 0.0;
2564	G4double pi = 3.141592653589793238e0;
2565
2566	// fundamental constants
2567	// hydrogen-atom mass excess
2568	mh = 7.289034;
2569
2570	// neutron mass excess
2571	mn = 8.071431;
2572
2573	// electronic charge squared
2574	esq = 1.4399764;
2575
2576	// constants from considerations other than nucl. masses
2577	// electronic binding
2578	ael = 1.433e-5;
2579
2580	// proton rms radius
2581	rp = 0.8;
2582
2583	// nuclear radius constant
2584	r0 = 1.16;
2585
2586	// range of yukawa-plus-expon. potential
2587	ay = 0.68;
2588
2589	// range of yukawa function used to generate
2590	// nuclear charge distribution
2591	aden= 0.70;
2592
2593	// constants from considering odd-even mass differences
2594	// average pairing gap
2595	rmac= 4.80;
2596
2597	// neutron-proton interaction
2598	h = 6.6;
2599
2600	// wigner constant
2601	w = 30.0;
2602
2603	// adjusted parameters
2604	// volume energy
2605	av = 16.00126;
2606
2607	// volume asymmetry
2608	kv = 1.92240;
2609
2610	// surface energy
2611	as = 21.18466;
2612
2613	// surface asymmetry
2614	ks = 2.345;
2615	// a^0 constant
2616	a0 = 2.615;
2617
2618	// charge asymmetry
2619	ca = 0.10289;
2620
2621	// we will account for deformation by using the microscopic
2622	// corrections tabulated from p. 68ff */
2623	bs = 1.0;
2624
2625	z = double(iz);
2626	a = double(ia);
2627	in = ia - iz;
2628	n = double(in);
2629	dn = rmac*bs/std::pow(n,(1.0/3.0));
2630	dp = rmac*bs/std::pow(z,(1.0/3.0));
2631	dpn = h/bs/std::pow(a,(2.0/3.0));
2632
2633	c1 = 3.0/5.0*esq/r0;
2634	c4 = 5.0/4.0std::pow((3.0/(2.0pi)),(2.0/3.0)) * c1;
2635	kf = std::pow((9.0piz/(4.0*a)),(1.0/3.0))/r0;
2636
2637	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));
2638	i = (n-z)/a;
2639
2640	x0 = r0 * std::pow(a,(1.0/3.0)) / ay;
2641	y0 = r0 * std::pow(a,(1.0/3.0)) / aden;
2642
2643	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);
2644
2645	b3 = 1.0 - 5.0/std::pow(y0,2) * (1.0 - 15.0/(8.0y0) + 21.0/(8.0 std::pow(y0,3))
2646	- 3.0/4.0 * (1.0 + 9.0/(2.0*y0) + 7.0/std::pow(y0,2)
2647	+ 7.0/(2.0 * std::pow(y0,3))) * std::exp(-2.0*y0));
2648
2649	// now calulation of total binding energy a.j. 16.7.96
2650
2651	efl = -1.0 * av(1.0 - kvii)a + as(1.0 - ksii)b1 * std::pow(a,(2.0/3.0)) + a0
2652	+ c1zzb3/std::pow(a,(1.0/3.0)) - c4std::pow(z,(4.0/3.0))/std::pow(a,(1.e0/3.e0))
2653	+ fstd::pow(z,2)/a -ca(n-z) - ael * std::pow(z,(2.39e0));
2654
2655	if ((in == iz) && (mod(in,2) == 1) && (mod(iz,2) == 1)) {
2656	// n and z odd and equal
2657	efl = efl + w*(utilabs(i)+1.e0/a);
2658	}
2659	else {
2660	efl= efl + w* utilabs(i);
2661	}
2662
2663	// pairing is made optional
2664	if (optshp >= 2) {
2665	// average pairing
2666	if ((mod(in,2) == 1) && (mod(iz,2) == 1)) {
2667	efl = efl - dpn;
2668	}
2669	if (mod(in,2) == 1) {
2670	efl = efl + dn;
2671	}
2672	if (mod(iz,2) == 1) {
2673	efl = efl + dp;
2674	}
2675	// end if for pairing term
2676	}
2677
2678	if (flag != 0) {
2679	eflmacResult = (0.5(dn + dp) - 0.5dpn);
2680	}
2681	else {
2682	eflmacResult = efl;
2683	}
2684
2685	return eflmacResult;
2686	}
2687
2688	void G4Abla::appariem(G4double a, G4double z, G4double *del)
2689	{
2690	// CALCUL DE LA CORRECTION, DUE A L'APPARIEMENT, DE L'ENERGIE DE
2691	// LIAISON D'UN NOYAU
2692	// PROCEDURE FOR CALCULATING THE PAIRING CORRECTION TO THE BINDING
2693	// ENERGY OF A SPECIFIC NUCLEUS
2694
2695	double para = 0.0, parz = 0.0;
2696	// A MASS NUMBER
2697	// Z NUCLEAR CHARGE
2698	// PARA HELP VARIABLE FOR PARITY OF A
2699	// PARZ HELP VARIABLE FOR PARITY OF Z
2700	// DEL PAIRING CORRECTION
2701
2702	parite(a, &para);
2703
2704	if (para < 0.0) {
2705	(*del) = 0.0;
2706	}
2707	else {
2708	parite(z, &parz);
2709	if (parz > 0.0) {
2710	(*del) = -12.0/std::sqrt(a);
2711	}
2712	else {
2713	(*del) = 12.0/std::sqrt(a);
2714	}
2715	}
2716	}
2717
2718	void G4Abla::parite(G4double n, G4double *par)
2719	{
2720	// CALCUL DE LA PARITE DU NOMBRE N
2721	//
2722	// PROCEDURE FOR CALCULATING THE PARITY OF THE NUMBER N.
2723	// RETURNS -1 IF N IS ODD AND +1 IF N IS EVEN
2724
2725	G4double n1 = 0.0, n2 = 0.0, n3 = 0.0;
2726
2727	// N NUMBER TO BE TESTED
2728	// N1,N2 HELP VARIABLES
2729	// PAR HELP VARIABLE FOR PARITY OF N
2730
2731	n3 = double(idnint(n));
2732	n1 = n3/2.0;
2733	n2 = n1 - dint(n1);
2734
2735	if (n2 > 0.0) {
2736	(*par) = -1.0;
2737	}
2738	else {
2739	(*par) = 1.0;
2740	}
2741	}
2742
2743	G4double G4Abla::tau(G4double bet, G4double homega, G4double ef, G4double t)
2744	{
2745	// INPUT : BET, HOMEGA, EF, T
2746	// OUTPUT: TAU - RISE TIME IN WHICH THE FISSION WIDTH HAS REACHED
2747	// 90 PERCENT OF ITS FINAL VALUE
2748	//
2749	// BETA - NUCLEAR VISCOSITY
2750	// HOMEGA - CURVATURE OF POTENTIAL
2751	// EF - FISSION BARRIER
2752	// T - NUCLEAR TEMPERATURE
2753
2754	G4double tauResult = 0.0;
2755
2756	G4double tlim = 8.e0 * ef;
2757	if (t > tlim) {
2758	t = tlim;
2759	}
2760
2761	// modified bj and khs 6.1.2000:
2762	if (bet/(std::sqrt(2.0)10.0(homega/6.582122)) <= 1.0) {
2763	tauResult = std::log(10.0ef/t)/(bet1.0e21);
2764	}
2765	else {
2766	tauResult = std::log(10.0ef/t)/ (2.0std::pow((10.0homega/6.582122),2))(bet*1.0e-21);
2767	} //end if
2768
2769	return tauResult;
2770	}
2771
2772	G4double G4Abla::cram(G4double bet, G4double homega)
2773	{
2774	// INPUT : BET, HOMEGA NUCLEAR VISCOSITY + CURVATURE OF POTENTIAL
2775	// OUTPUT: KRAMERS FAKTOR - REDUCTION OF THE FISSION PROBABILITY
2776	// INDEPENDENT OF EXCITATION ENERGY
2777
2778	G4double rel = bet/(20.0*homega/6.582122);
2779	G4double cramResult = std::sqrt(1.0 + std::pow(rel,2)) - rel;
2780	// limitation introduced 6.1.2000 by khs
2781
2782	if (cramResult > 1.0) {
2783	cramResult = 1.0;
2784	}
2785
2786	return cramResult;
2787	}
2788
2789	G4double G4Abla::bipol(int iflag, G4double y)
2790	{
2791	// CALCULATION OF THE SURFACE BS OR CURVATURE BK OF A NUCLEUS
2792	// RELATIVE TO THE SPHERICAL CONFIGURATION
2793	// BASED ON MYERS, DROPLET MODEL FOR ARBITRARY SHAPES
2794
2795	// INPUT: IFLAG - 0/1 BK/BS CALCULATION
2796	// Y - (1 - X) COMPLEMENT OF THE FISSILITY
2797
2798	// LINEAR INTERPOLATION OF BS BK TABLE
2799
2800	int i = 0;
2801
2802	G4double bipolResult = 0.0;
2803
2804	const int bsbkSize = 54;
2805
2806	G4double bk[bsbkSize] = {0.0, 1.00000,1.00087,1.00352,1.00799,1.01433,1.02265,1.03306,
2807	1.04576,1.06099,1.07910,1.10056,1.12603,1.15651,1.19348,
2808	1.23915,1.29590,1.35951,1.41013,1.44103,1.46026,1.47339,
2809	1.48308,1.49068,1.49692,1.50226,1.50694,1.51114,1.51502,
2810	1.51864,1.52208,1.52539,1.52861,1.53177,1.53490,1.53803,
2811	1.54117,1.54473,1.54762,1.55096,1.55440,1.55798,1.56173,
2812	1.56567,1.56980,1.57413,1.57860,1.58301,1.58688,1.58688,
2813	1.58688,1.58740,1.58740, 0.0}; //Zeroes at bk[0], and at the end added by PK
2814
2815	G4double bs[bsbkSize] = {0.0, 1.00000,1.00086,1.00338,1.00750,1.01319,
2816	1.02044,1.02927,1.03974,
2817	1.05195,1.06604,1.08224,1.10085,1.12229,1.14717,1.17623,1.20963,
2818	1.24296,1.26532,1.27619,1.28126,1.28362,1.28458,1.28477,1.28450,
2819	1.28394,1.28320,1.28235,1.28141,1.28042,1.27941,1.27837,1.27732,
2820	1.27627,1.27522,1.27418,1.27314,1.27210,1.27108,1.27006,1.26906,
2821	1.26806,1.26707,1.26610,1.26514,1.26418,1.26325,1.26233,1.26147,
2822	1.26147,1.26147,1.25992,1.25992, 0.0};
2823
2824	i = idint(y/(2.0e-02)) + 1;
2825
2826	if(i >= bsbkSize) {
2827	if(verboseLevel > 2) {
2828	G4cout <<"G4Abla error: index i = " << i << " is greater than array size permits." << G4endl;
2829	}
2830	bipolResult = 0.0;
2831	}
2832	else {
2833	if (iflag == 1) {
2834	bipolResult = bs[i] + (bs[i+1] - bs[i])/2.0e-02 * (y - 2.0e-02*(i - 1));
2835	}
2836	else {
2837	bipolResult = bk[i] + (bk[i+1] - bk[i])/2.0e-02 * (y - 2.0e-02*(i - 1));
2838	}
2839	}
2840
2841	return bipolResult;
2842	}
2843
2844	void G4Abla::barfit(G4int iz, G4int ia, G4int il, G4double sbfis, G4double segs, G4double *selmax)
2845	{
2846	// 2223 C VERSION FOR 32BIT COMPUTER
2847	// 2224 C THIS SUBROUTINE RETURNS THE BARRIER HEIGHT BFIS, THE
2848	// 2225 C GROUND-STATE ENERGY SEGS, IN MEV, AND THE ANGULAR MOMENTUM
2849	// 2226 C AT WHICH THE FISSION BARRIER DISAPPEARS, LMAX, IN UNITS OF
2850	// 2227 C H-BAR, WHEN CALLED WITH INTEGER AGUMENTS IZ, THE ATOMIC
2851	// 2228 C NUMBER, IA, THE ATOMIC MASS NUMBER, AND IL, THE ANGULAR
2852	// 2229 C MOMENTUM IN UNITS OF H-BAR. (PLANCK'S CONSTANT DIVIDED BY
2853	// 2230 C 2*PI).
2854	// 2231 C
2855	// 2232 C THE FISSION BARRIER FO IL = 0 IS CALCULATED FROM A 7TH
2856	// 2233 C ORDER FIT IN TWO VARIABLES TO 638 CALCULATED FISSION
2857	// 2234 C BARRIERS FOR Z VALUES FROM 20 TO 110. THESE 638 BARRIERS ARE
2858	// 2235 C FIT WITH AN RMS DEVIATION OF 0.10 MEV BY THIS 49-PARAMETER
2859	// 2236 C FUNCTION.
2860	// 2237 C IF BARFIT IS CALLED WITH (IZ,IA) VALUES OUTSIDE THE RANGE OF
2861	// 2238 C THE BARRIER HEIGHT IS SET TO 0.0, AND A MESSAGE IS PRINTED
2862	// 2239 C ON THE DEFAULT OUTPUT FILE.
2863	// 2240 C
2864	// 2241 C FOR IL VALUES NOT EQUAL TO ZERO, THE VALUES OF L AT WHICH
2865	// 2242 C THE BARRIER IS 80% AND 20% OF THE L=0 VALUE ARE RESPECTIVELY
2866	// 2243 C FIT TO 20-PARAMETER FUNCTIONS OF Z AND A, OVER A MORE
2867	// 2244 C RESTRICTED RANGE OF A VALUES, THAN IS THE CASE FOR L = 0.
2868	// 2245 C THE VALUE OF L WHERE THE BARRIER DISAPPEARS, LMAX IS FIT TO
2869	// 2246 C A 24-PARAMETER FUNCTION OF Z AND A, WITH THE SAME RANGE OF
2870	// 2247 C Z AND A VALUES AS L-80 AND L-20.
2871	// 2248 C ONCE AGAIN, IF AN (IZ,IA) PAIR IS OUTSIDE OF THE RANGE OF
2872	// 2249 C VALIDITY OF THE FIT, THE BARRIER VALUE IS SET TO 0.0 AND A
2873	// 2250 C MESSAGE IS PRINTED. THESE THREE VALUES (BFIS(L=0),L-80, AND
2874	// 2251 C L-20) AND THE CONSTRINTS OF BFIS = 0 AND D(BFIS)/DL = 0 AT
2875	// 2252 C L = LMAX AND L=0 LEAD TO A FIFTH-ORDER FIT TO BFIS(L) FOR
2876	// 2253 C L>L-20. THE FIRST THREE CONSTRAINTS LEAD TO A THIRD-ORDER FIT
2877	// 2254 C FOR THE REGION L < L-20.
2878	// 2255 C
2879	// 2256 C THE GROUND STATE ENERGIES ARE CALCULATED FROM A
2880	// 2257 C 120-PARAMETER FIT IN Z, A, AND L TO 214 GROUND-STATE ENERGIES
2881	// 2258 C FOR 36 DIFFERENT Z AND A VALUES.
2882	// 2259 C (THE RANGE OF Z AND A IS THE SAME AS FOR L-80, L-20, AND
2883	// 2260 C L-MAX)
2884	// 2261 C
2885	// 2262 C THE CALCULATED BARRIERS FROM WHICH THE FITS WERE MADE WERE
2886	// 2263 C CALCULATED IN 1983-1984 BY A. J. SIERK OF LOS ALAMOS
2887	// 2264 C NATIONAL LABORATORY GROUP T-9, USING YUKAWA-PLUS-EXPONENTIAL
2888	// 2265 C G4DOUBLE FOLDED NUCLEAR ENERGY, EXACT COULOMB DIFFUSENESS
2889	// 2266 C CORRECTIONS, AND DIFFUSE-MATTER MOMENTS OF INERTIA.
2890	// 2267 C THE PARAMETERS OF THE MODEL R-0 = 1.16 FM, AS 21.13 MEV,
2891	// 2268 C KAPPA-S = 2.3, A = 0.68 FM.
2892	// 2269 C THE DIFFUSENESS OF THE MATTER AND CHARGE DISTRIBUTIONS USED
2893	// 2270 C CORRESPONDS TO A SURFACE DIFFUSENESS PARAMETER (DEFINED BY
2894	// 2271 C MYERS) OF 0.99 FM. THE CALCULATED BARRIERS FOR L = 0 ARE
2895	// 2272 C ACCURATE TO A LITTLE LESS THAN 0.1 MEV; THE OUTPUT FROM
2896	// 2273 C THIS SUBROUTINE IS A LITTLE LESS ACCURATE. WORST ERRORS MAY BE
2897	// 2274 C AS LARGE AS 0.5 MEV; CHARACTERISTIC UNCERTAINY IS IN THE RANGE
2898	// 2275 C OF 0.1-0.2 MEV. THE RMS DEVIATION OF THE GROUND-STATE FIT
2899	// 2276 C FROM THE 214 INPUT VALUES IS 0.20 MEV. THE MAXIMUM ERROR
2900	// 2277 C OCCURS FOR LIGHT NUCLEI IN THE REGION WHERE THE GROUND STATE
2901	// 2278 C IS PROLATE, AND MAY BE GREATER THAN 1.0 MEV FOR VERY NEUTRON
2902	// 2279 C DEFICIENT NUCLEI, WITH L NEAR LMAX. FOR MOST NUCLEI LIKELY TO
2903	// 2280 C BE ENCOUNTERED IN REAL EXPERIMENTS, THE MAXIMUM ERROR IS
2904	// 2281 C CLOSER TO 0.5 MEV, AGAIN FOR LIGHT NUCLEI AND L NEAR LMAX.
2905	// 2282 C
2906	// 2283 C WRITTEN BY A. J. SIERK, LANL T-9
2907	// 2284 C VERSION 1.0 FEBRUARY, 1984
2908	// 2285 C
2909	// 2286 C THE FOLLOWING IS NECESSARY FOR 32-BIT MACHINES LIKE DEC VAX,
2910	// 2287 C IBM, ETC
2911
2912	G4double pa[7],pz[7],pl[10];
2913	for(G4int init_i = 0; init_i < 7; init_i++) {
2914	pa[init_i] = 0.0;
2915	pz[init_i] = 0.0;
2916	}
2917	for(G4int init_i = 0; init_i < 10; init_i++) {
2918	pl[init_i] = 0.0;
2919	}
2920
2921	G4double a = 0.0, z = 0.0, amin = 0.0, amax = 0.0, amin2 = 0.0;
2922	G4double amax2 = 0.0, aa = 0.0, zz = 0.0, bfis = 0.0;
2923	G4double bfis0 = 0.0, ell = 0.0, el = 0.0, egs = 0.0, el80 = 0.0, el20 = 0.0;
2924	G4double elmax = 0.0, sel80 = 0.0, sel20 = 0.0, x = 0.0, y = 0.0, q = 0.0, qa = 0.0, qb = 0.0;
2925	G4double aj = 0.0, ak = 0.0, a1 = 0.0, a2 = 0.0;
2926
2927	G4int i = 0, j = 0, k = 0, m = 0;
2928	G4int l = 0;
2929
2930	G4double emncof[4][5] = {{-9.01100e+2,-1.40818e+3, 2.77000e+3,-7.06695e+2, 8.89867e+2},
2931	{1.35355e+4,-2.03847e+4, 1.09384e+4,-4.86297e+3,-6.18603e+2},
2932	{-3.26367e+3, 1.62447e+3, 1.36856e+3, 1.31731e+3, 1.53372e+2},
2933	{7.48863e+3,-1.21581e+4, 5.50281e+3,-1.33630e+3, 5.05367e-2}};
2934
2935	G4double elmcof[4][5] = {{1.84542e+3,-5.64002e+3, 5.66730e+3,-3.15150e+3, 9.54160e+2},
2936	{-2.24577e+3, 8.56133e+3,-9.67348e+3, 5.81744e+3,-1.86997e+3},
2937	{2.79772e+3,-8.73073e+3, 9.19706e+3,-4.91900e+3, 1.37283e+3},
2938	{-3.01866e+1, 1.41161e+3,-2.85919e+3, 2.13016e+3,-6.49072e+2}};
2939
2940	G4double emxcof[4][6] = {{9.43596e4,-2.241997e5,2.223237e5,-1.324408e5,4.68922e4,-8.83568e3},
2941	{-1.655827e5,4.062365e5,-4.236128e5,2.66837e5,-9.93242e4,1.90644e4},
2942	{1.705447e5,-4.032e5,3.970312e5,-2.313704e5,7.81147e4,-1.322775e4},
2943	{-9.274555e4,2.278093e5,-2.422225e5,1.55431e5,-5.78742e4,9.97505e3}};
2944
2945	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},
2946	{-1.13269453e+6, 2.97764590e+6,-4.54326326e+6, 3.00464870e+6, -1.44989274e+6,-1.02026610e+5, 6.27959815e+4},
2947	{1.37543304e+6,-3.65808988e+6, 5.47798999e+6,-3.78109283e+6, 1.84131765e+6, 1.53669695e+4,-6.96817834e+4},
2948	{-8.56559835e+5, 2.48872266e+6,-4.07349128e+6, 3.12835899e+6, -1.62394090e+6, 1.19797378e+5, 4.25737058e+4},
2949	{3.28723311e+5,-1.09892175e+6, 2.03997269e+6,-1.77185718e+6, 9.96051545e+5,-1.53305699e+5,-1.12982954e+4},
2950	{4.15850238e+4, 7.29653408e+4,-4.93776346e+5, 6.01254680e+5, -4.01308292e+5, 9.65968391e+4,-3.49596027e+3},
2951	{-1.82751044e+5, 3.91386300e+5,-3.03639248e+5, 1.15782417e+5, -4.24399280e+3,-6.11477247e+3, 3.66982647e+2}};
2952
2953	const G4int sizex = 5;
2954	const G4int sizey = 6;
2955	const G4int sizez = 4;
2956
2957	G4double egscof[sizey][sizey][sizez];
2958
2959	G4double egs1[sizey][sizex] = {{1.927813e5, 7.666859e5, 6.628436e5, 1.586504e5,-7.786476e3},
2960	{-4.499687e5,-1.784644e6,-1.546968e6,-4.020658e5,-3.929522e3},
2961	{4.667741e5, 1.849838e6, 1.641313e6, 5.229787e5, 5.928137e4},
2962	{-3.017927e5,-1.206483e6,-1.124685e6,-4.478641e5,-8.682323e4},
2963	{1.226517e5, 5.015667e5, 5.032605e5, 2.404477e5, 5.603301e4},
2964	{-1.752824e4,-7.411621e4,-7.989019e4,-4.175486e4,-1.024194e4}};
2965
2966	G4double egs2[sizey][sizex] = {{-6.459162e5,-2.903581e6,-3.048551e6,-1.004411e6,-6.558220e4},
2967	{1.469853e6, 6.564615e6, 6.843078e6, 2.280839e6, 1.802023e5},
2968	{-1.435116e6,-6.322470e6,-6.531834e6,-2.298744e6,-2.639612e5},
2969	{8.665296e5, 3.769159e6, 3.899685e6, 1.520520e6, 2.498728e5},
2970	{-3.302885e5,-1.429313e6,-1.512075e6,-6.744828e5,-1.398771e5},
2971	{4.958167e4, 2.178202e5, 2.400617e5, 1.167815e5, 2.663901e4}};
2972
2973	G4double egs3[sizey][sizex] = {{3.117030e5, 1.195474e6, 9.036289e5, 6.876190e4,-6.814556e4},
2974	{-7.394913e5,-2.826468e6,-2.152757e6,-2.459553e5, 1.101414e5},
2975	{7.918994e5, 3.030439e6, 2.412611e6, 5.228065e5, 8.542465e3},
2976	{-5.421004e5,-2.102672e6,-1.813959e6,-6.251700e5,-1.184348e5},
2977	{2.370771e5, 9.459043e5, 9.026235e5, 4.116799e5, 1.001348e5},
2978	{-4.227664e4,-1.738756e5,-1.795906e5,-9.292141e4,-2.397528e4}};
2979
2980	G4double egs4[sizey][sizex] = {{-1.072763e5,-5.973532e5,-6.151814e5, 7.371898e4, 1.255490e5},
2981	{2.298769e5, 1.265001e6, 1.252798e6,-2.306276e5,-2.845824e5},
2982	{-2.093664e5,-1.100874e6,-1.009313e6, 2.705945e5, 2.506562e5},
2983	{1.274613e5, 6.190307e5, 5.262822e5,-1.336039e5,-1.115865e5},
2984	{-5.715764e4,-2.560989e5,-2.228781e5,-3.222789e3, 1.575670e4},
2985	{1.189447e4, 5.161815e4, 4.870290e4, 1.266808e4, 2.069603e3}};
2986
2987	for(i = 0; i < sizey; i++) {
2988	for(j = 0; j < sizex; j++) {
2989	// egscof[i][j][0] = egs1[i][j];
2990	// egscof[i][j][1] = egs2[i][j];
2991	// egscof[i][j][2] = egs3[i][j];
2992	// egscof[i][j][3] = egs4[i][j];
2993	egscof[i][j][0] = egs1[i][j];
2994	egscof[i][j][1] = egs2[i][j];
2995	egscof[i][j][2] = egs3[i][j];
2996	egscof[i][j][3] = egs4[i][j];
2997	}
2998	}
2999
3000	// the program starts here
3001	if (iz < 19 \|\| iz > 111) {
3002	goto barfit900;
3003	}
3004
3005	if(iz > 102 && il > 0) {
3006	goto barfit902;
3007	}
3008
3009	z=double(iz);
3010	a=double(ia);
3011	el=double(il);
3012	amin= 1.2e0z + 0.01e0z*z;
3013	amax= 5.8e0z - 0.024e0z*z;
3014
3015	if(a < amin \|\| a > amax) {
3016	goto barfit910;
3017	}
3018
3019	// angul.mom.zero barrier
3020	aa=2.5e-3*a;
3021	zz=1.0e-2*z;
3022	ell=1.0e-2*el;
3023	bfis0 = 0.0;
3024	lpoly(zz,7,pz);
3025	lpoly(aa,7,pa);
3026
3027	for(i = 0; i < 7; i++) { //do 10 i=1,7
3028	for(j = 0; j < 7; j++) { //do 10 j=1,7
3029	bfis0=bfis0+elzcof[j][i]pz[i]pa[j];
3030	//bfis0=bfis0+elzcof[i][j]pz[j]pa[i];
3031	}
3032	}
3033
3034	bfis=bfis0;
3035
3036	(*sbfis)=bfis;
3037	egs=0.0;
3038	(*segs)=egs;
3039
3040	// values of l at which the barrier
3041	// is 20%(el20) and 80%(el80) of l=0 value
3042	amin2 = 1.4e0z + 0.009e0z*z;
3043	amax2 = 20.e0 + 3.0e0*z;
3044
3045	if((a < amin2-5.e0 \|\| a > amax2+10.e0) && il > 0) {
3046	goto barfit920;
3047	}
3048
3049	lpoly(zz,5,pz);
3050	lpoly(aa,4,pa);
3051	el80=0.0;
3052	el20=0.0;
3053	elmax=0.0;
3054
3055	for(i = 0; i < 4; i++) {
3056	for(j = 0; j < 5; j++) {
3057	// el80 = el80 + elmcof[j][i]pz[j]pa[i];
3058	// el20 = el20 + emncof[j][i]pz[j]pa[i];
3059	el80 = el80 + elmcof[i][j]pz[j]pa[i];
3060	el20 = el20 + emncof[i][j]pz[j]pa[i];
3061	}
3062	}
3063
3064	sel80 = el80;
3065	sel20 = el20;
3066
3067	// value of l (elmax) where barrier disapp.
3068	lpoly(zz,6,pz);
3069	lpoly(ell,9,pl);
3070
3071	for(i = 0; i < 4; i++) { //do 30 i= 1,4
3072	for(j = 0; j < 6; j++) { //do 30 j=1,6
3073	//elmax = elmax + emxcof[j][i]pz[j]pa[i];
3074	// elmax = elmax + emxcof[j][i]pz[i]pa[j];
3075	elmax = elmax + emxcof[i][j]pz[j]pa[i];
3076	}
3077	}
3078
3079	(*selmax)=elmax;
3080
3081	// value of barrier at ang.mom. l
3082	if(il < 1){
3083	return;
3084	}
3085
3086	x = sel20/(*selmax);
3087	y = sel80/(*selmax);
3088
3089	if(el <= sel20) {
3090	// low l
3091	q = 0.2/(std::pow(sel20,2)std::pow(sel80,2)(sel20-sel80));
3092	qa = q(4.0std::pow(sel80,3) - std::pow(sel20,3));
3093	qb = -q(4.0std::pow(sel80,2) - std::pow(sel20,2));
3094	bfis = bfis(1.0 + qastd::pow(el,2) + qb*std::pow(el,3));
3095	}
3096	else {
3097	// high l
3098	aj = (-20.0std::pow(x,5) + 25.e0std::pow(x,4) - 4.0)std::pow((y-1.0),2)y*y;
3099	ak = (-20.0std::pow(y,5) + 25.0std::pow(y,4) - 1.0) * std::pow((x-1.0),2)xx;
3100	q = 0.2/(std::pow((y-x)((1.0-x)(1.0-y)xy),2));
3101	qa = q(ajy - ak*x);
3102	qb = -q(aj(2.0y + 1.0) - ak(2.0*x + 1.0));
3103	z = el/(*selmax);
3104	a1 = 4.0std::pow(z,5) - 5.0std::pow(z,4) + 1.0;
3105	a2 = qa(2.e0z + 1.e0);
3106	bfis=bfis(a1 + (z - 1.e0)(a2 + qbz)zz(z - 1.e0));
3107	}
3108
3109	if(bfis <= 0.0) {
3110	bfis=0.0;
3111	}
3112
3113	if(el > (*selmax)) {
3114	bfis=0.0;
3115	}
3116	(*sbfis)=bfis;
3117
3118	// now calculate rotating ground state energy
3119	if(el > (*selmax)) {
3120	return;
3121	}
3122
3123	for(k = 0; k < 4; k++) {
3124	for(l = 0; l < 6; l++) {
3125	for(m = 0; m < 5; m++) {
3126	//egs = egs + egscof[l][m][k]pz[l]pa[k]pl[2m-1];
3127	egs = egs + egscof[l][m][k]pz[l]pa[k]pl[2m];
3128	// egs = egs + egscof[m][l][k]pz[l]pa[k]pl[2m-1];
3129	}
3130	}
3131	}
3132
3133	(*segs)=egs;
3134	if((*segs) < 0.0) {
3135	(*segs)=0.0;
3136	}
3137
3138	return;
3139
3140	barfit900: //continue
3141	(*sbfis)=0.0;
3142	// for z<19 sbfis set to 1.0e3
3143	if (iz < 19) {
3144	(*sbfis) = 1.0e3;
3145	}
3146	(*segs)=0.0;
3147	(*selmax)=0.0;
3148	return;
3149
3150	barfit902:
3151	(*sbfis)=0.0;
3152	(*segs)=0.0;
3153	(*selmax)=0.0;
3154	return;
3155
3156	barfit910:
3157	(*sbfis)=0.0;
3158	(*segs)=0.0;
3159	(*selmax)=0.0;
3160	return;
3161
3162	barfit920:
3163	(*sbfis)=0.0;
3164	(*segs)=0.0;
3165	(*selmax)=0.0;
3166	return;
3167	}
3168
3169	G4double G4Abla::expohaz(G4int k, G4double T)
3170	{
3171	// TIRAGE ALEATOIRE DANS UNE EXPONENTIELLLE : Y=EXP(-X/T)
3172
3173	return (-1.0Tstd::log(haz(k)));
3174	}
3175
3176	G4double G4Abla::fd(G4double E)
3177	{
3178	// DISTRIBUTION DE MAXWELL
3179
3180	return (E*std::exp(-E));
3181	}
3182
3183	G4double G4Abla::f(G4double E)
3184	{
3185	// FONCTION INTEGRALE DE FD(E)
3186	return (1.0 - (E + 1.0) * std::exp(-E));
3187	}
3188
3189	G4double G4Abla::fmaxhaz(G4double T)
3190	{
3191	// tirage aleatoire dans une maxwellienne
3192	// t : temperature
3193	//
3194	// declaration des variables
3195	//
3196
3197	const int pSize = 101;
3198	static G4double p[pSize];
3199
3200	// ial generateur pour le cascade (et les iy pour eviter les correlations)
3201	static G4int i = 0;
3202	static G4int itest = 0;
3203	// programme principal
3204
3205	// calcul des p(i) par approximation de newton
3206	p[pSize-1] = 8.0;
3207	G4double x = 0.1;
3208	G4double x1 = 0.0;
3209	G4double y = 0.0;
3210
3211	if (itest == 1) {
3212	goto fmaxhaz120;
3213	}
3214
3215	for(i = 1; i <= 99; i++) {
3216	fmaxhaz20:
3217	x1 = x - (f(x) - double(i)/100.0)/fd(x);
3218	x = x1;
3219	if (std::fabs(f(x) - double(i)/100.0) < 1e-5) {
3220	goto fmaxhaz100;
3221	}
3222	goto fmaxhaz20;
3223	fmaxhaz100:
3224	p[i] = x;
3225	} //end do
3226
3227	// itest = 1;
3228	itest = 0;
3229	// tirage aleatoire et calcul du x correspondant
3230	// par regression lineaire
3231	fmaxhaz120:
3232	y = randomGenerator->getRandom();
3233	// standardRandom(&y, &(hazard->igraine[17]));
3234	i = nint(y*100);
3235
3236	// 2590 c ici on evite froidement les depassements de tableaux....(a.b. 3/9/99)
3237	if(i == 0) {
3238	goto fmaxhaz120;
3239	}
3240
3241	if (i == 1) {
3242	x = p[i]y100;
3243	}
3244	else {
3245	x = (p[i] - p[i-1])(y100 - i) + p[i];
3246	}
3247
3248	return(x*T);
3249	}
3250
3251	G4double G4Abla::pace2(G4double a, G4double z)
3252	{
3253	// PACE2
3254	// Cette fonction retourne le defaut de masse du noyau A,Z en MeV
3255	// RÃ©visÃ©e pour a, z flottants 25/4/2002 =
3256
3257	G4double pace2 = 0.0;
3258
3259	G4int ii = idint(a+0.5);
3260	G4int jj = idint(z+0.5);
3261
3262	if(ii <= 0 \|\| jj < 0) {
3263	pace2=0.;
3264	return pace2;
3265	}
3266
3267	if(jj > 300) {
3268	pace2=0.0;
3269	}
3270	else {
3271	pace2=pace->dm[ii][jj];
3272	}
3273	pace2=pace2/1000.;
3274
3275	if(pace->dm[ii][jj] == 0.) {
3276	if(ii < 12) {
3277	pace2=-500.;
3278	}
3279	else {
3280	guet(&a, &z, &pace2);
3281	pace2=pace2-ii*931.5;
3282	pace2=pace2/1000.;
3283	}
3284	}
3285
3286	return pace2;
3287	}
3288
3289	void G4Abla::guet(G4double x_par, G4double z_par, G4double *find_par)
3290	{
3291	// TABLE DE MASSES ET FORMULE DE MASSE TIRE DU PAPIER DE BRACK-GUET
3292	// Gives the theoritical value for mass excess...
3293	// RÃ©visÃ©e pour x, z flottants 25/4/2002
3294
3295	//real*8 x,z
3296	// dimension q(0:50,0:70)
3297	G4double x = (*x_par);
3298	G4double z = (*z_par);
3299	G4double find = (*find_par);
3300
3301	const G4int qrows = 50;
3302	const G4int qcols = 70;
3303	G4double q[qrows][qcols];
3304	for(G4int init_i = 0; init_i < qrows; init_i++) {
3305	for(G4int init_j = 0; init_j < qcols; init_j++) {
3306	q[init_i][init_j] = 0.0;
3307	}
3308	}
3309
3310	G4int ix=G4int(std::floor(x+0.5));
3311	G4int iz=G4int(std::floor(z+0.5));
3312	G4double zz = iz;
3313	G4double xx = ix;
3314	find = 0.0;
3315	G4double avol = 15.776;
3316	G4double asur = -17.22;
3317	G4double ac = -10.24;
3318	G4double azer = 8.0;
3319	G4double xjj = -30.03;
3320	G4double qq = -35.4;
3321	G4double c1 = -0.737;
3322	G4double c2 = 1.28;
3323
3324	if(ix <= 7) {
3325	q[0][1]=939.50;
3326	q[1][1]=938.21;
3327	q[1][2]=1876.1;
3328	q[1][3]=2809.39;
3329	q[2][4]=3728.34;
3330	q[2][3]=2809.4;
3331	q[2][5]=4668.8;
3332	q[2][6]=5606.5;
3333	q[3][5]=4669.1;
3334	q[3][6]=5602.9;
3335	q[3][7]=6535.27;
3336	q[4][6]=5607.3;
3337	q[4][7]=6536.1;
3338	q[5][7]=6548.3;
3339	find=q[iz][ix];
3340	}
3341	else {
3342	G4double xneu=xx-zz;
3343	G4double si=(xneu-zz)/xx;
3344	G4double x13=std::pow(xx,.333);
3345	G4double ee1=c1zzzz/x13;
3346	G4double ee2=c2zzzz/xx;
3347	G4double aux=1.+(9.*xjj/4./qq/x13);
3348	G4double ee3=xjjxxsi*si/aux;
3349	G4double ee4=avolxx+asur(std::pow(xx,.666))+ac*x13+azer;
3350	G4double tota = ee1 + ee2 + ee3 + ee4;
3351	find = 939.55xneu+938.77zz - tota;
3352	}
3353
3354	(*x_par) = x;
3355	(*z_par) = z;
3356	(*find_par) = find;
3357	}
3358
3359	// SUBROUTINE TRANSLAB(GAMREM,ETREM,CSREM,NOPART,NDEC)
3360	void G4Abla::translab(G4double gamrem, G4double etrem, G4double csrem[4], G4int nopart, G4int ndec)
3361	{
3362	// c Ce subroutine transforme dans un repere 1 les impulsions pcv des
3363	// c particules acv, zcv et de cosinus directeurs xcv, ycv, zcv calculees
3364	// c dans un repere 2.
3365	// c La transformation de lorentz est definie par GAMREM (gamma) et
3366	// c ETREM (eta). La direction du repere 2 dans 1 est donnees par les
3367	// c cosinus directeurs ALREM,BEREM,GAREM (axe oz du repere 2).
3368	// c L'axe oy(2) est fixe par le produit vectoriel oz(1)*oz(2).
3369	// c Le calcul est fait pour les particules de NDEC a iv du common volant.
3370	// C Resultats dans le NTUPLE (common VAR_NTP) decale de NOPART (cascade).
3371
3372	// REAL*8 GAMREM,ETREM,ER,PLABI(3),PLABF(3),R(3,3)
3373	// real*8 MASSE,PTRAV2,CSREM(3),UMA,MELEC,EL
3374	// real*4 acv,zpcv,pcv,xcv,ycv,zcv
3375	// common/volant/acv(200),zpcv(200),pcv(200),xcv(200),
3376	// s ycv(200),zcv(200),iv
3377
3378	// parameter (max=250)
3379	// real*4 EXINI,ENERJ,BIMPACT,PLAB,TETLAB,PHILAB,ESTFIS
3380	// integer AVV,ZVV,JREMN,KFIS,IZFIS,IAFIS
3381	// common/VAR_NTP/MASSINI,MZINI,EXINI,MULNCASC,MULNEVAP,
3382	// +MULNTOT,BIMPACT,JREMN,KFIS,ESTFIS,IZFIS,IAFIS,NTRACK,
3383	// +ITYPCASC(max),AVV(max),ZVV(max),ENERJ(max),PLAB(max),
3384	// +TETLAB(max),PHILAB(max)
3385
3386	// DATA UMA,MELEC/931.4942,0.511/
3387
3388	// C Matrice de rotation dans le labo:
3389	G4double avv = 0.0, zvv = 0.0, enerj = 0.0, plab = 0.0, tetlab = 0.0, philab = 0.0;
3390	G4int itypcasc = 0;
3391	G4double sitet = std::sqrt(std::pow(csrem[1],2)+std::pow(csrem[2],2));
3392	G4double cstet = 0.0, siphi = 0.0, csphi = 0.0;
3393	G4double R[4][4];
3394	for(G4int init_i = 0; init_i < 4; init_i++) {
3395	for(G4int init_j = 0; init_j < 4; init_j++) {
3396	R[init_i][init_j] = 0.0;
3397	}
3398	}
3399	if(verboseLevel > 1)
3400	G4cout <<"csrem = " << csrem[1] << ", " << csrem[2] << ", " << csrem[3] << ")" << G4endl;
3401	if(sitet > 1.0e-6) { //then
3402	cstet = csrem[3];
3403	siphi = csrem[2]/sitet;
3404	csphi = csrem[1]/sitet;
3405
3406	R[1][1] = cstet*csphi;
3407	R[1][2] = -siphi;
3408	R[1][3] = sitet*csphi;
3409	R[2][1] = cstet*siphi;
3410	R[2][2] = csphi;
3411	R[2][3] = sitet*siphi;
3412	R[3][1] = -sitet;
3413	R[3][2] = 0.0;
3414	R[3][3] = cstet;
3415	}
3416	else {
3417	R[1][1] = 1.0;
3418	R[1][2] = 0.0;
3419	R[1][3] = 0.0;
3420	R[2][1] = 0.0;
3421	R[2][2] = 1.0;
3422	R[2][3] = 0.0;
3423	R[3][1] = 0.0;
3424	R[3][2] = 0.0;
3425	R[3][3] = 1.0;
3426	} //endif
3427
3428	G4int intp = 0;
3429	G4double el = 0.0;
3430	G4double masse = 0.0;
3431	G4double er = 0.0;
3432	G4double plabi[4];
3433	G4double ptrav2 = 0.0;
3434	G4double plabf[4];
3435	G4double bidon = 0.0;
3436	for(G4int init_i = 0; init_i < 4; init_i++) {
3437	plabi[init_i] = 0.0;
3438	plabf[init_i] = 0.0;
3439	}
3440	ndec = 1;
3441	for(G4int i = ndec; i <= volant->iv; i++) { //do i=ndec,iv
3442	intp = i + nopart;
3443	if(volant->copied[i]) continue; // Avoid double copying
3444	#ifdef USE_LEGACY_CODE
3445	varntp->ntrack = varntp->ntrack + 1;
3446	#endif
3447	if(nint(volant->acv[i]) == 0 && nint(volant->zpcv[i]) == 0) {
3448	if(verboseLevel > -1) {
3449	G4cout << __FILE__ << ":" << __LINE__ << " Error: Particles with A = 0 Z = 0 detected! " << G4endl;
3450	}
3451	continue;
3452	}
3453	if(varntp->ntrack >= VARNTPSIZE) {
3454	if(verboseLevel > 2) {
3455	G4cout <<"Error! Output data structure not big enough!" << G4endl;
3456	}
3457	}
3458	#ifdef USE_LEGACY_CODE
3459	varntp->avv[intp] = nint(volant->acv[i]);
3460	varntp->zvv[intp] = nint(volant->zpcv[i]);
3461	varntp->itypcasc[intp] = 0;
3462	#else
3463	avv = nint(volant->acv[i]);
3464	zvv = nint(volant->zpcv[i]);
3465	itypcasc = 0;
3466	#endif
3467	// transformation de lorentz remnan --> labo:
3468	#ifdef USE_LEGACY_CODE
3469	if (varntp->avv[intp] == -1) { //then
3470	#else
3471	if(avv == -1) {
3472	#endif
3473	masse = 138.00; // cugnon
3474	// c if (avv(intp).eq.1) masse=938.2796 !cugnon
3475	// c if (avv(intp).eq.4) masse=3727.42 !ok
3476	}
3477	else {
3478	mglms(double(volant->acv[i]),double(volant->zpcv[i]),0, &el);
3479	masse = volant->zpcv[i]938.27 + (volant->acv[i] - volant->zpcv[i])939.56 + el;
3480	} //end if
3481
3482	er = std::sqrt(std::pow(volant->pcv[i],2) + std::pow(masse,2));
3483	plabi[1] = volant->pcv[i]*(volant->xcv[i]);
3484	plabi[2] = volant->pcv[i]*(volant->ycv[i]);
3485	plabi[3] = eretrem + gamrem(volant->pcv[i])*(volant->zcv[i]);
3486
3487	ptrav2 = std::pow(plabi[1],2) + std::pow(plabi[2],2) + std::pow(plabi[3],2);
3488	#ifdef USE_LEGACY_CODE
3489	varntp->plab[intp] = std::sqrt(ptrav2);
3490	#else
3491	plab = std::sqrt(ptrav2);
3492	#endif
3493	#ifdef USE_LEGACY_CODE
3494	if(std::abs(varntp->plab[intp] - 122.009) < 1.0e-3) {
3495	G4cout <<__FILE__ << ":" << __LINE__ << " Error: varntp->plab["<< intp <<"] = " << varntp->plab[intp] << G4endl;
3496
3497	volant->dump();
3498	varntp->dump();
3499	#
3500	G4cout <<"varntp->plab[intp] = " << varntp->plab[intp] << G4endl;
3501	G4cout <<"ndec (starting index for loop) = " << ndec << G4endl;
3502	G4cout <<"volant->iv (stopping index for loop) = " << volant->iv << G4endl;
3503	G4cout <<"i (current position) = " << i << G4endl;
3504	G4cout <<"intp (index for writing to varntp) = " << intp << G4endl;
3505	// exit(0);
3506	}
3507	#endif
3508
3509	#ifdef USE_LEGACY_CODE
3510	varntp->enerj[intp] = std::sqrt(ptrav2 + std::pow(masse,2)) - masse;
3511	#else
3512	enerj = std::sqrt(ptrav2 + std::pow(masse,2)) - masse;
3513	#endif
3514	// Rotation dans le labo:
3515	for(G4int j = 1; j <= 3; j++) { //do j=1,3
3516	plabf[j] = 0.0;
3517	for(G4int k = 1; k <= 3; k++) { //do k=1,3
3518	//plabf[j] = plabf[j] + R[k][j]*plabi[k]; // :::Fixme::: (indices?)
3519	plabf[j] = plabf[j] + R[j][k]*plabi[k]; // :::Fixme::: (indices?)
3520	} // end do
3521	} // end do
3522	// C impulsions dans le nouveau systeme copiees dans /volant/
3523	#ifdef USE_LEGACY_CODE
3524	volant->pcv[i] = varntp->plab[intp];
3525	#else
3526	volant->pcv[i] = plab;
3527	#endif
3528	ptrav2 = std::sqrt(std::pow(plabf[1],2) + std::pow(plabf[2],2) + std::pow(plabf[3],2));
3529	if(ptrav2 >= 1.0e-6) { //then
3530	volant->xcv[i] = plabf[1]/ptrav2;
3531	volant->ycv[i] = plabf[2]/ptrav2;
3532	volant->zcv[i] = plabf[3]/ptrav2;
3533	}
3534	else {
3535	volant->xcv[i] = 1.0;
3536	volant->ycv[i] = 0.0;
3537	volant->zcv[i] = 0.0;
3538	} //endif
3539	// impulsions dans le nouveau systeme copiees dans /VAR_NTP/
3540	#ifdef USE_LEGACY_CODE
3541	if(varntp->plab[intp] >= 1.0e-6) { //then
3542	#else
3543	if(plab >= 1.0e-6) { //then
3544	#endif
3545	#ifdef USE_LEGACY_CODE
3546	bidon = plabf[3]/(varntp->plab[intp]);
3547	#else
3548	bidon = plabf[3]/plab;
3549	#endif
3550	if(bidon > 1.0) {
3551	bidon = 1.0;
3552	}
3553	if(bidon < -1.0) {
3554	bidon = -1.0;
3555	}
3556	#ifdef USE_LEGACY_CODE
3557	varntp->tetlab[intp] = std::acos(bidon);
3558	sitet = std::sin(varntp->tetlab[intp]);
3559	varntp->philab[intp] = std::atan2(plabf[2],plabf[1]);
3560	varntp->tetlab[intp] = varntp->tetlab[intp]*57.2957795;
3561	varntp->philab[intp] = varntp->philab[intp]*57.2957795;
3562	#else
3563	tetlab = std::acos(bidon);
3564	sitet = std::sin(tetlab);
3565	philab = std::atan2(plabf[2],plabf[1]);
3566	tetlab = tetlab*57.2957795;
3567	philab = philab*57.2957795;
3568	#endif
3569	}
3570	else {
3571	#ifdef USE_LEGACY_CODE
3572	varntp->tetlab[intp] = 90.0;
3573	varntp->philab[intp] = 0.0;
3574	#else
3575	tetlab = 90.0;
3576	philab = 0.0;
3577	#endif
3578	} // endif
3579	volant->copied[i] = true;
3580	#ifndef USE_LEGACY_CODE
3581	varntp->addParticle(avv, zvv, enerj, plab, tetlab, philab);
3582	#else
3583	G4cout <<__FILE__ << ":" << __LINE__ << " volant -> varntp: " << " intp = " << intp << "Avv = " << varntp->avv[intp] << " Zvv = " << varntp->zvv[intp] << " Plab = " << varntp->plab[intp] << G4endl;
3584	G4cout <<__FILE__ << ":" << __LINE__ << " volant index i = " << i << " varnpt index intp = " << intp << G4endl;
3585	#endif
3586	} // end do
3587	}
3588	// C-------------------------------------------------------------------------
3589
3590	// SUBROUTINE TRANSLABPF(MASSE1,T1,P1,CTET1,PHI1,GAMREM,ETREM,R,
3591	// s PLAB1,GAM1,ETA1,CSDIR)
3592	void G4Abla::translabpf(G4double masse1, G4double t1, G4double p1, G4double ctet1,
3593	G4double phi1, G4double gamrem, G4double etrem, G4double R[][4],
3594	G4double plab1, G4double gam1, G4double *eta1, G4double csdir[])
3595	{
3596	// C Calcul de l'impulsion du PF (PLAB1, cos directeurs CSDIR(3)) dans le
3597	// C systeme remnant et des coefs de Lorentz GAM1,ETA1 de passage
3598	// c du systeme PF --> systeme remnant.
3599	// c
3600	// C Input: MASSE1, T1 (energie cinetique), CTET1,PHI1 (cosTHETA et PHI)
3601	// C (le PF dans le systeme du Noyau de Fission (NF)).
3602	// C GAMREM,ETREM les coefs de Lorentz systeme NF --> syst remnant,
3603	// C R(3,3) la matrice de rotation systeme NF--> systeme remnant.
3604	// C
3605	// C
3606	// REAL*8 MASSE1,T1,P1,CTET1,PHI1,GAMREM,ETREM,R(3,3),
3607	// s PLAB1,GAM1,ETA1,CSDIR(3),ER,SITET,PLABI(3),PLABF(3)
3608
3609	G4double er = t1 + masse1;
3610
3611	G4double sitet = std::sqrt(1.0 - std::pow(ctet1,2));
3612
3613	G4double plabi[4];
3614	G4double plabf[4];
3615	for(G4int init_i = 0; init_i < 4; init_i++) {
3616	plabi[init_i] = 0.0;
3617	plabf[init_i] = 0.0;
3618	}
3619
3620	// C ----Transformation de Lorentz Noyau fissionnant --> Remnant:
3621	plabi[1] = p1sitetstd::cos(phi1);
3622	plabi[2] = p1sitetstd::sin(phi1);
3623	plabi[3] = eretrem + gamremp1*ctet1;
3624
3625	// C ----Rotation du syst Noyaut Fissionant vers syst remnant:
3626	for(G4int j = 1; j <= 3; j++) { // do j=1,3
3627	plabf[j] = 0.0;
3628	for(G4int k = 1; k <= 3; k++) { //do k=1,3
3629	plabf[j] = plabf[j] + R[j][k]*plabi[k];
3630	//plabf[j] = plabf[j] + R[k][j]*plabi[k];
3631	} //end do
3632	} //end do
3633	// C ----Cosinus directeurs et coefs de la transf de Lorentz dans le
3634	// c nouveau systeme:
3635	(*plab1) = std::pow(plabf[1],2) + std::pow(plabf[2],2) + std::pow(plabf[3],2);
3636	(gam1) = std::sqrt(std::pow(masse1,2) + (plab1))/masse1;
3637	(plab1) = std::sqrt((plab1));
3638	(eta1) = (plab1)/masse1;
3639
3640	if((*plab1) <= 1.0e-6) { //then
3641	csdir[1] = 0.0;
3642	csdir[2] = 0.0;
3643	csdir[3] = 1.0;
3644	}
3645	else {
3646	for(G4int i = 1; i <= 3; i++) { //do i=1,3
3647	csdir[i] = plabf[i]/(*plab1);
3648	} // end do
3649	} //endif
3650	}
3651
3652	// SUBROUTINE LOR_AB(GAM,ETA,Ein,Pin,Eout,Pout)
3653	void G4Abla::lorab(G4double gam, G4double eta, G4double ein, G4double pin[],
3654	G4double *eout, G4double pout[])
3655	{
3656	// C Transformation de lorentz brute pour vÃ©rifs.
3657	// C P(3) = P_longitudinal (transformÃ©)
3658	// C P(1) et P(2) = P_transvers (non transformÃ©s)
3659	// DIMENSION Pin(3),Pout(3)
3660	// REAL*8 GAM,ETA,Ein
3661
3662	pout[1] = pin[1];
3663	pout[2] = pin[2];
3664	(eout) = gamein + eta*pin[3];
3665	pout[3] = etaein + gampin[3];
3666	}
3667
3668	// SUBROUTINE ROT_AB(R,Pin,Pout)
3669	void G4Abla::rotab(G4double R[4][4], G4double pin[4], G4double pout[4])
3670	{
3671	// C Rotation d'un vecteur
3672	// DIMENSION Pin(3),Pout(3)
3673	// REAL*8 R(3,3)
3674
3675	for(G4int i = 1; i <= 3; i++) { // do i=1,3
3676	pout[i] = 0.0;
3677	for(G4int j = 1; j <= 3; j++) { //do j=1,3
3678	pout[i] = pout[i] + R[i][j]*pin[j];
3679	//pout[i] = pout[i] + R[j][i]*pin[j];
3680	} // enddo
3681	} //enddo
3682	}
3683
3684	// Methods related to the internal ABLA random number generator. In
3685	// the future the random number generation must be factored into its
3686	// own class
3687
3688	void G4Abla::standardRandom(G4double rndm, G4long seed)
3689	{
3690	(seed) = (seed); // Avoid warning during compilation.
3691	// Use Geant4 G4UniformRand
3692	(*rndm) = randomGenerator->getRandom();
3693	}
3694
3695	G4double G4Abla::haz(G4int k)
3696	{
3697	const G4int pSize = 110;
3698	static G4double p[pSize];
3699	static G4long ix = 0, i = 0;
3700	static G4double x = 0.0, y = 0.0, a = 0.0, haz = 0.0;
3701	// k =< -1 on initialise
3702	// k = -1 c'est reproductible
3703	// k < -1 \|\| k > -1 ce n'est pas reproductible
3704
3705	// Zero is invalid random seed. Set proper value from our random seed collection:
3706	if(ix == 0) {
3707	ix = hazard->ial;
3708	}
3709
3710	if (k <= -1) { //then
3711	if(k == -1) { //then
3712	ix = 0;
3713	}
3714	else {
3715	x = 0.0;
3716	y = secnds(int(x));
3717	ix = int(y * 100 + 43543000);
3718	if(mod(ix,2) == 0) {
3719	ix = ix + 1;
3720	}
3721	}
3722
3723	// Here we are using random number generator copied from INCL code
3724	// instead of the CERNLIB one! This causes difficulties for
3725	// automatic testing since the random number generators, and thus
3726	// the behavior of the routines in C++ and FORTRAN versions is no
3727	// longer exactly the same!
3728	x = randomGenerator->getRandom();
3729	// standardRandom(&x, &ix);
3730	for(G4int i = 0; i < pSize; i++) { //do i=1,110
3731	p[i] = randomGenerator->getRandom();
3732	// standardRandom(&(p[i]), &ix);
3733	}
3734	a = randomGenerator->getRandom();
3735	standardRandom(&a, &ix);
3736	k = 0;
3737	}
3738
3739	i = nint(100*a)+1;
3740	haz = p[i];
3741	a = randomGenerator->getRandom();
3742	// standardRandom(&a, &ix);
3743	p[i] = a;
3744
3745	hazard->ial = ix;
3746	return haz;
3747	}
3748
3749	// Utilities
3750
3751	G4double G4Abla::min(G4double a, G4double b)
3752	{
3753	if(a < b) {
3754	return a;
3755	}
3756	else {
3757	return b;
3758	}
3759	}
3760
3761	G4int G4Abla::min(G4int a, G4int b)
3762	{
3763	if(a < b) {
3764	return a;
3765	}
3766	else {
3767	return b;
3768	}
3769	}
3770
3771	G4double G4Abla::max(G4double a, G4double b)
3772	{
3773	if(a > b) {
3774	return a;
3775	}
3776	else {
3777	return b;
3778	}
3779	}
3780
3781	G4int G4Abla::max(G4int a, G4int b)
3782	{
3783	if(a > b) {
3784	return a;
3785	}
3786	else {
3787	return b;
3788	}
3789	}
3790
3791	G4int G4Abla::nint(G4double number)
3792	{
3793	G4double intpart = 0.0;
3794	G4double fractpart = 0.0;
3795	fractpart = std::modf(number, &intpart);
3796	if(number == 0) {
3797	return 0;
3798	}
3799	if(number > 0) {
3800	if(fractpart < 0.5) {
3801	return int(std::floor(number));
3802	}
3803	else {
3804	return int(std::ceil(number));
3805	}
3806	}
3807	if(number < 0) {
3808	if(fractpart < -0.5) {
3809	return int(std::floor(number));
3810	}
3811	else {
3812	return int(std::ceil(number));
3813	}
3814	}
3815
3816	return int(std::floor(number));
3817	}
3818
3819	G4int G4Abla::secnds(G4int x)
3820	{
3821	time_t mytime;
3822	tm *mylocaltime;
3823
3824	time(&mytime);
3825	mylocaltime = localtime(&mytime);
3826
3827	if(x == 0) {
3828	return(mylocaltime->tm_hour6060 + mylocaltime->tm_min*60 + mylocaltime->tm_sec);
3829	}
3830	else {
3831	return(mytime - x);
3832	}
3833	}
3834
3835	G4int G4Abla::mod(G4int a, G4int b)
3836	{
3837	if(b != 0) {
3838	return (a - (a/b)*b);
3839	}
3840	else {
3841	return 0;
3842	}
3843	}
3844
3845	G4double G4Abla::dmod(G4double a, G4double b)
3846	{
3847	if(b != 0) {
3848	return (a - (a/b)*b);
3849	}
3850	else {
3851	return 0.0;
3852	}
3853	}
3854
3855	G4double G4Abla::dint(G4double a)
3856	{
3857	G4double value = 0.0;
3858
3859	if(a < 0.0) {
3860	value = double(std::ceil(a));
3861	}
3862	else {
3863	value = double(std::floor(a));
3864	}
3865
3866	return value;
3867	}
3868
3869	G4int G4Abla::idint(G4double a)
3870	{
3871	G4int value = 0;
3872
3873	if(a < 0) {
3874	value = int(std::ceil(a));
3875	}
3876	else {
3877	value = int(std::floor(a));
3878	}
3879
3880	return value;
3881	}
3882
3883	G4int G4Abla::idnint(G4double value)
3884	{
3885	G4double valueCeil = int(std::ceil(value));
3886	G4double valueFloor = int(std::floor(value));
3887
3888	if(std::fabs(value - valueCeil) < std::fabs(value - valueFloor)) {
3889	return int(valueCeil);
3890	}
3891	else {
3892	return int(valueFloor);
3893	}
3894	}
3895
3896	G4double G4Abla::dmin1(G4double a, G4double b, G4double c)
3897	{
3898	if(a < b && a < c) {
3899	return a;
3900	}
3901	if(b < a && b < c) {
3902	return b;
3903	}
3904	if(c < a && c < b) {
3905	return c;
3906	}
3907	return a;
3908	}
3909
3910	G4double G4Abla::utilabs(G4double a)
3911	{
3912	if(a > 0) {
3913	return a;
3914	}
3915	if(a < 0) {
3916	return (-1*a);
3917	}
3918	if(a == 0) {
3919	return a;
3920	}
3921
3922	return a;
3923	}

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

Download in other formats: