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

Last change on this file since 1201 was 1196, checked in by garnier, 16 years ago

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