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

Last change on this file since 1344 was 1340, checked in by garnier, 15 years ago

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