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

Last change on this file since 880 was 819, checked in by garnier, 17 years ago

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