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