Context Navigation

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

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

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

Download in other formats: