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