Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

G4Abla.cc

Last change on this file was 1347, checked in by garnier, 14 years ago
geant4 tag 9.4
File size: 140.2 KB

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

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

Context Navigation

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

Download in other formats: