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

Last change on this file since 1036 was 962, checked in by garnier, 17 years ago

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