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