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

Last change on this file since 1357 was 1347, checked in by garnier, 15 years ago

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