00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033 #include <time.h>
00034
00035 #include "G4Abla.hh"
00036 #include "G4InclAblaDataFile.hh"
00037 #include "Randomize.hh"
00038 #include <assert.h>
00039
00040 G4Abla::G4Abla()
00041 {
00042 ilast = 0;
00043 }
00044
00045 G4Abla::G4Abla(G4Hazard *hazard, G4Volant *volant)
00046 {
00047 verboseLevel = 0;
00048 ilast = 0;
00049 volant = volant;
00050 volant->iv = 0;
00051 hazard = hazard;
00052
00053 varntp = new G4VarNtp();
00054 pace = new G4Pace();
00055 ald = new G4Ald();
00056 ablamain = new G4Ablamain();
00057 emdpar = new G4Emdpar();
00058 eenuc = new G4Eenuc();
00059 ec2sub = new G4Ec2sub();
00060 ecld = new G4Ecld();
00061 fb = new G4Fb();
00062 fiss = new G4Fiss();
00063 opt = new G4Opt();
00064 }
00065
00066 G4Abla::G4Abla(G4Hazard *aHazard, G4Volant *aVolant, G4VarNtp *aVarntp)
00067 {
00068 verboseLevel = 0;
00069 ilast = 0;
00070 volant = aVolant;
00071 volant->iv = 0;
00072 hazard = aHazard;
00073 varntp = aVarntp;
00074 varntp->ntrack = 0;
00075
00076 pace = new G4Pace();
00077 ald = new G4Ald();
00078 ablamain = new G4Ablamain();
00079 emdpar = new G4Emdpar();
00080 eenuc = new G4Eenuc();
00081 ec2sub = new G4Ec2sub();
00082 ecld = new G4Ecld();
00083 fb = new G4Fb();
00084 fiss = new G4Fiss();
00085 opt = new G4Opt();
00086 }
00087
00088 G4Abla::~G4Abla()
00089 {
00090 delete pace;
00091 delete ald;
00092 delete ablamain;
00093 delete emdpar;
00094 delete eenuc;
00095 delete ec2sub;
00096 delete ecld;
00097 delete fb;
00098 delete fiss;
00099 delete opt;
00100 }
00101
00102
00103
00104
00105
00106
00107
00108
00109 void G4Abla::breakItUp(G4double nucleusA, G4double nucleusZ, G4double nucleusMass, G4double excitationEnergy,
00110 G4double angularMomentum, G4double recoilEnergy, G4double momX, G4double momY, G4double momZ,
00111 G4int eventnumber)
00112 {
00113 const G4double uma = 931.4942;
00114 const G4double melec = 0.511;
00115 const G4double fmp = 938.27231;
00116 const G4double fmn = 939.56563;
00117
00118 G4double alrem = 0.0, berem = 0.0, garem = 0.0;
00119 G4double R[4][4];
00120 G4double csdir1[4];
00121 G4double csdir2[4];
00122 G4double csrem[4];
00123 G4double pfis_rem[4];
00124 G4double pf1_rem[4];
00125 for(G4int init_i = 0; init_i < 4; init_i++) {
00126 csdir1[init_i] = 0.0;
00127 csdir2[init_i] = 0.0;
00128 csrem[init_i] = 0.0;
00129 pfis_rem[init_i] = 0.0;
00130 pf1_rem[init_i] = 0.0;
00131 for(G4int init_j = 0; init_j < 4; init_j++) {
00132 R[init_i][init_j] = 0.0;
00133 }
00134 }
00135
00136 G4double plab1 = 0.0, gam1 = 0.0, eta1 = 0.0;
00137 G4double plab2 = 0.0, gam2 = 0.0, eta2 = 0.0;
00138
00139 G4double sitet = 0.0;
00140 G4double stet1 = 0.0;
00141 G4double stet2 = 0.0;
00142
00143 G4int nbpevap = 0;
00144 G4int mempaw = 0, memiv = 0;
00145
00146 G4double e_evapo = 0.0;
00147 G4double el = 0.0;
00148 G4double fmcv = 0.0;
00149
00150 G4double aff1 = 0.0;
00151 G4double zff1 = 0.0;
00152 G4double eff1 = 0.0;
00153 G4double aff2 = 0.0;
00154 G4double zff2 = 0.0;
00155 G4double eff2 = 0.0;
00156
00157 G4double v1 = 0.0, v2 = 0.0;
00158
00159 G4double t2 = 0.0;
00160 G4double ctet1 = 0.0;
00161 G4double ctet2 = 0.0;
00162 G4double phi1 = 0.0;
00163 G4double phi2 = 0.0;
00164 G4double p2 = 0.0;
00165 G4double epf2_out = 0.0 ;
00166 G4int lma_pf1 = 0, lmi_pf1 = 0;
00167 G4int lma_pf2 = 0, lmi_pf2 = 0;
00168 G4int nopart = 0;
00169
00170 G4double cst = 0.0, sst = 0.0, csf = 0.0, ssf = 0.0;
00171
00172 G4double zf = 0.0, af = 0.0, mtota = 0.0, pleva = 0.0, pxeva = 0.0, pyeva = 0.0;
00173 G4int ff = 0;
00174 G4int inum = eventnumber;
00175 G4int inttype = 0;
00176 G4double esrem = excitationEnergy;
00177
00178 G4double aprf = nucleusA;
00179 G4double zprf = nucleusZ;
00180 G4double mcorem = nucleusMass;
00181 G4double ee = excitationEnergy;
00182 G4double jprf = angularMomentum;
00183
00184 G4double erecrem = recoilEnergy;
00185 G4double trem = 0.0;
00186 G4double pxrem = momX;
00187 G4double pyrem = momY;
00188 G4double pzrem = momZ;
00189
00190 G4double remmass = 0.0;
00191
00192 varntp->ntrack = 0;
00193
00194 volant->iv = 1;
00195
00196 G4double pcorem = std::sqrt(erecrem*(erecrem +2.*938.2796*nucleusA));
00197
00198
00199 if(esrem >= 1.0e-3) {
00200 evapora(zprf,aprf,ee,jprf, &zf, &af, &mtota, &pleva, &pxeva, &pyeva, &ff, &inttype, &inum);
00201 }
00202 else {
00203 ff = 0;
00204 zf = zprf;
00205 af = aprf;
00206 pxeva = pxrem;
00207 pyeva = pyrem;
00208 pleva = pzrem;
00209 }
00210
00211
00212
00213
00214 if (ff == 1) {
00215
00216
00217
00218
00219
00220 trem = double(erecrem);
00221 remmass = pace2(aprf,zprf) + aprf*uma - zprf*melec;
00222
00223
00224 varntp->kfis = 1;
00225 G4double gamrem = (remmass + trem)/remmass;
00226 G4double etrem = std::sqrt(trem*(trem + 2.0*remmass))/remmass;
00227
00228
00229
00230
00231
00232 remmass = pace2(aprf,zprf) + aprf*uma - zprf*melec+double(esrem);
00233
00234 el = 0.0;
00235 mglms(aprf,zprf,0,&el);
00236 remmass = zprf*fmp + (aprf-zprf)*fmn + el + double(esrem);
00237
00238 gamrem = std::sqrt(std::pow(pcorem,2) + std::pow(remmass,2))/remmass;
00239
00240 etrem = pcorem/remmass;
00241
00242
00243 alrem = pxrem/pcorem;
00244
00245 berem = pyrem/pcorem;
00246
00247 garem = pzrem/pcorem;
00248
00249
00250 csrem[0] = 0.0;
00251 csrem[1] = alrem;
00252 csrem[2] = berem;
00253 csrem[3] = garem;
00254
00255
00256 G4double bil_e = 0.0;
00257 G4double bil_px = 0.0;
00258 G4double bil_py = 0.0;
00259 G4double bil_pz = 0.0;
00260 G4double masse = 0.0;
00261
00262 for(G4int iloc = 1; iloc <= volant->iv; iloc++) {
00263
00264
00265
00266 mglms(double(volant->acv[iloc]),double(volant->zpcv[iloc]),0,&el);
00267
00268 masse = volant->zpcv[iloc]*fmp + (volant->acv[iloc] - volant->zpcv[iloc])*fmn + el;
00269
00270 bil_e = bil_e + std::sqrt(std::pow(volant->pcv[iloc],2) + std::pow(masse,2));
00271
00272 bil_px = bil_px + volant->pcv[iloc]*(volant->xcv[iloc]);
00273 bil_py = bil_py + volant->pcv[iloc]*(volant->ycv[iloc]);
00274 bil_pz = bil_pz + volant->pcv[iloc]*(volant->zcv[iloc]);
00275 }
00276
00277
00278 G4int ndec = 1;
00279
00280 if(volant->iv != 0) {
00281 if(verboseLevel > 2) {
00282 G4cout <<"varntp->ntrack = " << varntp->ntrack << G4endl;
00283 G4cout <<"1st Translab: Adding indices from " << ndec << " to " << volant->iv << G4endl;
00284 }
00285 nopart = varntp->ntrack - 1;
00286 translab(gamrem,etrem,csrem,nopart,ndec);
00287 if(verboseLevel > 2) {
00288 G4cout <<"Translab complete!" << G4endl;
00289 G4cout <<"varntp->ntrack = " << varntp->ntrack << G4endl;
00290 }
00291 }
00292 nbpevap = volant->iv;
00293
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303
00304
00305 fissionDistri(af,zf,ee,aff1,zff1,eff1,v1,aff2,zff2,eff2,v2);
00306
00307
00308 G4int na_f = int(std::floor(af + 0.5));
00309 G4int nz_f = int(std::floor(zf + 0.5));
00310 varntp->izfis = nz_f;
00311 varntp->iafis = na_f;
00312 G4int na_pf1 = int(std::floor(aff1 + 0.5));
00313 G4int nz_pf1 = int(std::floor(zff1 + 0.5));
00314 G4int na_pf2 = int(std::floor(aff2 + 0.5));
00315 G4int nz_pf2 = int(std::floor(zff2 + 0.5));
00316
00317 if((na_f != (na_pf1+na_pf2)) || (nz_f != (nz_pf1+nz_pf2))) {
00318 if(verboseLevel > 2) {
00319 G4cout <<"problemes arrondis dans la fission " << G4endl;
00320 G4cout << "af,zf,aff1,zff1,aff2,zff2" << G4endl;
00321 G4cout << af <<" , " << zf <<" , " << aff1 <<" , " << zff1 <<" , " << aff2 <<" , " << zff2 << G4endl;
00322 G4cout << "a,z,a1,z1,a2,z2 integer" << G4endl;
00323 G4cout << na_f <<" , " << nz_f <<" , " << na_pf1 <<" , " << nz_pf1 <<" , " << na_pf2 <<" , " << nz_pf2 << G4endl;
00324 }
00325 }
00326
00327
00328 G4int kboud = idnint(zf);
00329 G4int jboud = idnint(af-zf);
00330
00331 G4double ef = fb->efa[jboud][kboud];
00332
00333 varntp->estfis = ee + ef;
00334
00335
00336
00337
00338
00339
00340
00341
00342
00343
00344 mglms(af,zf,0,&el);
00345
00346 G4double massef = zf*fmp + (af - zf)*fmn + el + ee + ef;
00347
00348 mglms(double(aff1),double(zff1),0,&el);
00349
00350 G4double masse1 = zff1*fmp + (aff1-zff1)*fmn + el + eff1;
00351
00352 mglms(aff2,zff2,0,&el);
00353
00354 G4double masse2 = zff2*fmp + (aff2 - zff2)*fmn + el + eff2;
00355
00356
00357 G4double b = massef - masse1 - masse2;
00358 if(b < 0.0) {
00359 b=0.0;
00360 if(verboseLevel > 2) {
00361 G4cout <<"anomalie dans la fission: " << G4endl;
00362 G4cout << inum<< " , " << af<< " , " <<zf<< " , " <<massef<< " , " <<aff1<< " , " <<zff1<< " , " <<masse1<< " , " <<aff2<< " , " <<zff2<< " , " << masse2 << G4endl;
00363 }
00364 }
00365 G4double t1 = b*(b + 2.0*masse2)/(2.0*massef);
00366
00367 G4double p1 = std::sqrt(t1*(t1 + 2.0*masse1));
00368
00369
00370 G4double rndm;
00371 standardRandom(&rndm, &(hazard->igraine[13]));
00372 ctet1 = 2.0*rndm - 1.0;
00373 standardRandom(&rndm,&(hazard->igraine[9]));
00374 phi1 = rndm*2.0*3.141592654;
00375
00376
00377 G4double peva = std::pow(pxeva,2) + std::pow(pyeva,2) + std::pow(pleva,2);
00378 G4double gamfis = std::sqrt(std::pow(massef,2) + peva)/massef;
00379
00380 peva = std::sqrt(peva);
00381
00382 G4double etfis = peva/massef;
00383
00384 G4double epf1_in = 0.0;
00385 G4double epf1_out = 0.0;
00386
00387
00388 if(peva >= 1.0e-4) {
00389 sitet = std::sqrt(std::pow(pxeva,2)+std::pow(pyeva,2))/peva;
00390
00391 }
00392 if(sitet > 1.0e-5) {
00393 G4double cstet = pleva/peva;
00394 G4double siphi = pyeva/(sitet*peva);
00395 G4double csphi = pxeva/(sitet*peva);
00396
00397 R[1][1] = cstet*csphi;
00398 R[1][2] = -siphi;
00399 R[1][3] = sitet*csphi;
00400 R[2][1] = cstet*siphi;
00401 R[2][2] = csphi;
00402 R[2][3] = sitet*siphi;
00403 R[3][1] = -sitet;
00404 R[3][2] = 0.0;
00405 R[3][3] = cstet;
00406 }
00407 else {
00408 R[1][1] = 1.0;
00409 R[1][2] = 0.0;
00410 R[1][3] = 0.0;
00411 R[2][1] = 0.0;
00412 R[2][2] = 1.0;
00413 R[2][3] = 0.0;
00414 R[3][1] = 0.0;
00415 R[3][2] = 0.0;
00416 R[3][3] = 1.0;
00417 }
00418
00419
00420 if((zff1 <= 0.0) || (aff1 <= 0.0) || (aff1 < zff1)) {
00421 if(verboseLevel > 2) {
00422 G4cout <<"zf = " << zf <<" af = " << af <<"ee = " << ee <<"zff1 = " << zff1 <<"aff1 = " << aff1 << G4endl;
00423 }
00424 }
00425 else {
00426
00427 epf1_in = double(eff1);
00428 epf1_out = epf1_in;
00429
00430
00431
00432
00433 G4double zf1, af1, malpha1, ffpleva1, ffpxeva1, ffpyeva1;
00434 G4int ff1, ftype1;
00435 evapora(zff1, aff1, epf1_out, 0.0, &zf1, &af1, &malpha1, &ffpleva1,
00436 &ffpxeva1, &ffpyeva1, &ff1, &ftype1, &inum);
00437
00438
00439 volant->iv = volant->iv + 1;
00440
00441
00442 volant->acv[volant->iv] = af1;
00443 volant->zpcv[volant->iv] = zf1;
00444 if(verboseLevel > 2) {
00445 G4cout <<"Added fission fragment: a = " << volant->acv[volant->iv] << " z = " << volant->zpcv[volant->iv] << G4endl;
00446 }
00447 peva = std::sqrt(std::pow(ffpxeva1,2) + std::pow(ffpyeva1,2) + std::pow(ffpleva1,2));
00448
00449 volant->pcv[volant->iv] = peva;
00450 if(peva > 0.001) {
00451 volant->xcv[volant->iv] = ffpxeva1/peva;
00452 volant->ycv[volant->iv] = ffpyeva1/peva;
00453 volant->zcv[volant->iv] = ffpleva1/peva;
00454 }
00455 else {
00456 volant->xcv[volant->iv] = 1.0;
00457 volant->ycv[volant->iv] = 0.0;
00458 volant->zcv[volant->iv] = 0.0;
00459 }
00460
00461
00462 G4double bil1_e = 0.0;
00463 G4double bil1_px = 0.0;
00464 G4double bil1_py=0.0;
00465 G4double bil1_pz=0.0;
00466 for(G4int iloc = nbpevap + 1; iloc <= volant->iv; iloc++) {
00467
00468 mglms(volant->acv[iloc], volant->zpcv[iloc],0,&el);
00469 masse = volant->zpcv[iloc]*fmp + (volant->acv[iloc] - volant->zpcv[iloc])*fmn + el;
00470
00471 bil1_e = bil1_e + std::sqrt(std::pow(volant->pcv[iloc],2) + std::pow(masse,2));
00472
00473 bil1_px = bil1_px + volant->pcv[iloc]*(volant->xcv[iloc]);
00474 bil1_py = bil1_py + volant->pcv[iloc]*(volant->ycv[iloc]);
00475 bil1_pz = bil1_pz + volant->pcv[iloc]*(volant->zcv[iloc]);
00476 }
00477
00478
00479
00480 translabpf(masse1,t1,p1,ctet1,phi1,gamfis,etfis,R,&plab1,&gam1,&eta1,csdir1);
00481
00482
00483 if(verboseLevel > 2) {
00484 G4cout <<"2nd Translab (pf1 evap): Adding indices from " << nbpevap+1 << " to " << volant->iv << G4endl;
00485 }
00486 nopart = varntp->ntrack - 1;
00487 translab(gam1,eta1,csdir1,nopart,nbpevap+1);
00488 if(verboseLevel > 2) {
00489 G4cout <<"After translab call... varntp->ntrack = " << varntp->ntrack << G4endl;
00490 }
00491 memiv = nbpevap + 1;
00492 mempaw = nopart;
00493 lmi_pf1 = nopart + nbpevap + 1;
00494 lma_pf1 = nopart + volant->iv;
00495 nbpevap = volant->iv;
00496 }
00497
00498
00499
00500 if((zff2 <= 0.0) || (aff2 <= 0.0) || (aff2 <= zff2)) {
00501 if(verboseLevel > 2) {
00502 G4cout << zf << " " << af << " " << ee << " " << zff2 << " " << aff2 << G4endl;
00503 }
00504 }
00505 else {
00506
00507 G4double epf2_in = double(eff2);
00508 G4double epf2_out = epf2_in;
00509
00510
00511
00512
00513 G4double zf2, af2, malpha2, ffpleva2, ffpxeva2, ffpyeva2;
00514 G4int ff2, ftype2;
00515 evapora(zff2,aff2,epf2_out,0.0,&zf2,&af2,&malpha2,&ffpleva2,
00516 &ffpxeva2,&ffpyeva2,&ff2,&ftype2,&inum);
00517
00518 volant->iv = volant->iv + 1;
00519 volant->acv[volant->iv] = af2;
00520 volant->zpcv[volant->iv] = zf2;
00521 if(verboseLevel > 2) {
00522 G4cout <<"Added fission fragment: a = " << volant->acv[volant->iv] << " z = " << volant->zpcv[volant->iv] << G4endl;
00523 }
00524 peva = std::sqrt(std::pow(ffpxeva2,2) + std::pow(ffpyeva2,2) + std::pow(ffpleva2,2));
00525
00526 volant->pcv[volant->iv] = peva;
00527
00528 if(peva > 0.001) {
00529 volant->xcv[volant->iv] = ffpxeva2/peva;
00530 volant->ycv[volant->iv] = ffpyeva2/peva;
00531 volant->zcv[volant->iv] = ffpleva2/peva;
00532 }
00533 else {
00534 volant->xcv[volant->iv] = 1.0;
00535 volant->ycv[volant->iv] = 0.0;
00536 volant->zcv[volant->iv] = 0.0;
00537 }
00538
00539 G4double bil2_e = 0.0;
00540 G4double bil2_px = 0.0;
00541 G4double bil2_py = 0.0;
00542 G4double bil2_pz = 0.0;
00543
00544 for(G4int iloc = nbpevap + 1; iloc <= volant->iv; iloc++) {
00545 mglms(volant->acv[iloc],volant->zpcv[iloc],0,&el);
00546 masse = volant->zpcv[iloc]*fmp + (volant->acv[iloc] - volant->zpcv[iloc])*fmn + el;
00547 bil2_e = bil2_e + std::sqrt(std::pow(volant->pcv[iloc],2) + std::pow(masse,2));
00548
00549 bil2_px = bil2_px + volant->pcv[iloc]*(volant->xcv[iloc]);
00550 bil2_py = bil2_py + volant->pcv[iloc]*(volant->ycv[iloc]);
00551 bil2_pz = bil2_pz + volant->pcv[iloc]*(volant->zcv[iloc]);
00552 }
00553
00554
00555
00556 G4double t2 = b - t1;
00557
00558 ctet2 = -1.0*ctet1;
00559 assert(std::fabs(ctet2) <= 1.0);
00560
00561 phi2 = dmod(phi1+3.141592654,6.283185308);
00562
00563 G4double p2 = std::sqrt(t2*(t2+2.0*masse2));
00564
00565
00566
00567
00568
00569 translabpf(masse2,t2,p2,ctet2,phi2,gamfis,etfis,R,&plab2,&gam2,&eta2,csdir2);
00570
00571
00572
00573 if(verboseLevel > 2) {
00574 G4cout <<"3rd Translab (pf2 evap): Adding indices from " << nbpevap+1 << " to " << volant->iv << G4endl;
00575 }
00576 nopart = varntp->ntrack - 1;
00577 translab(gam2,eta2,csdir2,nopart,nbpevap+1);
00578 lmi_pf2 = nopart + nbpevap + 1;
00579 lma_pf2 = nopart + volant->iv;
00580 }
00581
00582
00583
00584
00585 for(G4int iloc = 1; iloc <= 3; iloc++) {
00586 pfis_rem[iloc] = 0.0;
00587 }
00588 G4double efis_rem, pfis_trav[4];
00589 lorab(gamfis,etfis,massef,pfis_rem,&efis_rem,pfis_trav);
00590 rotab(R,pfis_trav,pfis_rem);
00591
00592 stet1 = std::sqrt(1.0 - std::pow(ctet1,2));
00593
00594 pf1_rem[1] = p1*stet1*std::cos(phi1);
00595 pf1_rem[2] = p1*stet1*std::sin(phi1);
00596 pf1_rem[3] = p1*ctet1;
00597 G4double e1_rem;
00598 lorab(gamfis,etfis,masse1+t1,pf1_rem,&e1_rem,pfis_trav);
00599 rotab(R,pfis_trav,pf1_rem);
00600
00601 stet2 = std::sqrt(1.0 - std::pow(ctet2,2));
00602 assert(std::pow(ctet2,2) >= 0.0);
00603 assert(std::pow(ctet2,2) <= 1.0);
00604
00605
00606 G4double pf2_rem[4];
00607 G4double e2_rem;
00608 pf2_rem[1] = p2*stet2*std::cos(phi2);
00609 pf2_rem[2] = p2*stet2*std::sin(phi2);
00610 pf2_rem[3] = p2*ctet2;
00611 lorab(gamfis,etfis,masse2+t2,pf2_rem,&e2_rem,pfis_trav);
00612 rotab(R,pfis_trav,pf2_rem);
00613
00614 bil_e = remmass - efis_rem - bil_e;
00615 bil_px = bil_px + pfis_rem[1];
00616 bil_py = bil_py + pfis_rem[2];
00617 bil_pz = bil_pz + pfis_rem[3];
00618
00619
00620
00621
00622
00623
00624
00625 if((lma_pf1-lmi_pf1) != 0) {
00626 G4double bil_e_pf1 = e1_rem - epf1_out;
00627 G4double bil_px_pf1 = pf1_rem[1];
00628 G4double bil_py_pf1 = pf1_rem[2];
00629 G4double bil_pz_pf1 = pf1_rem[3];
00630 for(G4int ipf1 = lmi_pf1; ipf1 <= lma_pf1; ipf1++) {
00631 bil_e_pf1 = bil_e_pf1 - (std::pow(varntp->plab[ipf1],2) + std::pow(varntp->enerj[ipf1],2))/(2.0*(varntp->enerj[ipf1]));
00632 cst = std::cos(varntp->tetlab[ipf1]/57.2957795);
00633 sst = std::sin(varntp->tetlab[ipf1]/57.2957795);
00634 csf = std::cos(varntp->philab[ipf1]/57.2957795);
00635 ssf = std::sin(varntp->philab[ipf1]/57.2957795);
00636 bil_px_pf1 = bil_px_pf1 - varntp->plab[ipf1]*sst*csf;
00637 bil_py_pf1 = bil_py_pf1 - varntp->plab[ipf1]*sst*ssf;
00638 bil_pz_pf1 = bil_pz_pf1 - varntp->plab[ipf1]*cst;
00639 }
00640 }
00641
00642 if((lma_pf2-lmi_pf2) != 0) {
00643 G4double bil_e_pf2 = e2_rem - epf2_out;
00644 G4double bil_px_pf2 = pf2_rem[1];
00645 G4double bil_py_pf2 = pf2_rem[2];
00646 G4double bil_pz_pf2 = pf2_rem[3];
00647 for(G4int ipf2 = lmi_pf2; ipf2 <= lma_pf2; ipf2++) {
00648 bil_e_pf2 = bil_e_pf2 - (std::pow(varntp->plab[ipf2],2) + std::pow(varntp->enerj[ipf2],2))/(2.0*(varntp->enerj[ipf2]));
00649 G4double cst = std::cos(varntp->tetlab[ipf2]/57.2957795);
00650 G4double sst = std::sin(varntp->tetlab[ipf2]/57.2957795);
00651 G4double csf = std::cos(varntp->philab[ipf2]/57.2957795);
00652 G4double ssf = std::sin(varntp->philab[ipf2]/57.2957795);
00653 bil_px_pf2 = bil_px_pf2 - varntp->plab[ipf2]*sst*csf;
00654 bil_py_pf2 = bil_py_pf2 - varntp->plab[ipf2]*sst*ssf;
00655 bil_pz_pf2 = bil_pz_pf2 - varntp->plab[ipf2]*cst;
00656 }
00657 }
00658
00659
00660
00661
00662 if(verboseLevel > 2) {
00663 G4cout <<"4th Translab: Adding indices from " << memiv << " to " << volant->iv << G4endl;
00664 }
00665 translab(gamrem,etrem,csrem,mempaw,memiv);
00666
00667 }
00668 else {
00669
00670
00671
00672 varntp->kfis = 0;
00673 if(verboseLevel > 2) {
00674 G4cout <<"Evaporation without fission" << G4endl;
00675 }
00676 volant->iv = volant->iv + 1;
00677 volant->acv[volant->iv] = af;
00678 volant->zpcv[volant->iv] = zf;
00679 G4double peva = std::sqrt(std::pow(pxeva,2)+std::pow(pyeva,2)+std::pow(pleva,2));
00680
00681 volant->pcv[volant->iv] = peva;
00682 if(peva > 0.001) {
00683 volant->xcv[volant->iv] = pxeva/peva;
00684 volant->ycv[volant->iv] = pyeva/peva;
00685 volant->zcv[volant->iv] = pleva/peva;
00686 }
00687 else {
00688 volant->xcv[volant->iv] = 1.0;
00689 volant->ycv[volant->iv] = 0.0;
00690 volant->zcv[volant->iv] = 0.0;
00691 }
00692
00693
00694
00695
00696 trem = double(erecrem);
00697
00698
00699 remmass = mcorem;
00700
00701
00702 csrem[0] = 0.0;
00703 csrem[1] = alrem;
00704 csrem[2] = berem;
00705 csrem[3] = garem;
00706
00707
00708 for(G4int j = 1; j <= volant->iv; j++) {
00709 if(volant->acv[j] == 0) {
00710 if(verboseLevel > 2) {
00711 G4cout <<"volant->acv[" << j << "] = 0" << G4endl;
00712 G4cout <<"volant->iv = " << volant->iv << G4endl;
00713 }
00714 }
00715 if(volant->acv[j] > 0) {
00716 assert(volant->acv[j] != 0);
00717
00718 mglms(volant->acv[j],volant->zpcv[j],0,&el);
00719 fmcv = volant->zpcv[j]*fmp + (volant->acv[j] - volant->zpcv[j])*fmn + el;
00720 e_evapo = e_evapo + std::sqrt(std::pow(volant->pcv[j],2) + std::pow(fmcv,2));
00721
00722 }
00723 }
00724
00725
00726
00727
00728 remmass = e_evapo;
00729
00730 G4double gamrem = std::sqrt(std::pow(pcorem,2)+std::pow(remmass,2))/remmass;
00731
00732 G4double etrem = pcorem/remmass;
00733
00734 if(verboseLevel > 2) {
00735 G4cout <<"5th Translab (no fission): Adding indices from " << 1 << " to " << volant->iv << G4endl;
00736 }
00737 nopart = varntp->ntrack - 1;
00738 translab(gamrem,etrem,csrem,nopart,1);
00739
00740
00741 }
00742
00743 }
00744
00745
00746 void G4Abla::initEvapora()
00747 {
00748
00749
00750
00751
00752
00753
00754
00755
00756
00757
00758
00759
00760
00761
00762
00763
00764
00765
00766
00767
00768
00769
00770
00771
00772
00773
00774
00775
00776
00777
00778
00779
00780
00781
00782
00783
00784
00785
00786
00787
00788
00789
00790
00791
00792
00793
00794
00795
00796
00797
00798
00799
00800
00801
00802
00803
00804
00805
00806
00807
00808
00809
00810
00811
00812
00813
00814
00815
00816
00817
00818
00819
00820
00821
00822
00823
00824
00825
00826
00827
00828
00829
00830
00831
00832
00833
00834
00835
00836
00837
00838
00839
00840
00841
00842
00843
00844
00845
00846
00847
00848
00849
00850
00851
00852
00853
00854
00855
00856
00857
00858
00859
00860
00861
00862
00863
00864
00865
00866
00867
00868
00869
00870
00871
00872
00873
00874
00875
00876
00877
00878
00879
00880
00881
00882
00883
00884
00885
00886
00887
00888
00889
00890 fiss->ifis = 1;
00891
00892
00893 fiss->optshp = 0;
00894
00895
00896 opt->optemd = 1;
00897
00898 opt->optcha = 1;
00899
00900
00901 fiss->akap = 10.0;
00902
00903
00904 fiss->bet = 1.5;
00905
00906
00907 fiss->homega = 1.0;
00908
00909
00910 fiss->koeff = 1.;
00911
00912
00913 fiss->optcol = 0;
00914
00915
00916 fiss->optles = 0;
00917
00918 opt->eefac = 2.;
00919
00920 ald->optafan = 0;
00921
00922 ald->av = 0.073e0;
00923 ald->as = 0.095e0;
00924 ald->ak = 0.0e0;
00925
00926 if(verboseLevel > 3) {
00927 G4cout <<"ifis " << fiss->ifis << G4endl;
00928 G4cout <<"optshp " << fiss->optshp << G4endl;
00929 G4cout <<"optemd " << opt->optemd << G4endl;
00930 G4cout <<"optcha " << opt->optcha << G4endl;
00931 G4cout <<"akap " << fiss->akap << G4endl;
00932 G4cout <<"bet " << fiss->bet << G4endl;
00933 G4cout <<"homega " << fiss->homega << G4endl;
00934 G4cout <<"koeff " << fiss->koeff << G4endl;
00935 G4cout <<"optcol " << fiss->optcol << G4endl;
00936 G4cout <<"optles " << fiss->optles << G4endl;
00937 G4cout <<"eefac " << opt->eefac << G4endl;
00938 G4cout <<"optafan " << ald->optafan << G4endl;
00939 G4cout <<"av " << ald->av << G4endl;
00940 G4cout <<"as " << ald->as << G4endl;
00941 G4cout <<"ak " << ald->ak << G4endl;
00942 }
00943 fiss->optxfis = 1;
00944
00945 G4InclAblaDataFile *dataInterface = new G4InclAblaDataFile();
00946 if(dataInterface->readData() == true) {
00947 if(verboseLevel > 0) {
00948 G4cout <<"G4Abla: Datafiles read successfully." << G4endl;
00949 }
00950 }
00951 else {
00952 G4Exception("ERROR: Failed to read datafiles.");
00953 }
00954
00955 for(int z = 0; z < 98; z++) {
00956 for(int n = 0; n < 154; n++) {
00957 ecld->ecfnz[n][z] = 0.e0;
00958 ec2sub->ecnz[n][z] = dataInterface->getEcnz(n,z);
00959 ecld->ecgnz[n][z] = dataInterface->getEcnz(n,z);
00960 ecld->alpha[n][z] = dataInterface->getAlpha(n,z);
00961 ecld->vgsld[n][z] = dataInterface->getVgsld(n,z);
00962 }
00963 }
00964
00965 for(int z = 0; z < 500; z++) {
00966 for(int a = 0; a < 500; a++) {
00967 pace->dm[z][a] = dataInterface->getPace2(z,a);
00968 }
00969 }
00970
00971 delete dataInterface;
00972 }
00973
00974 void G4Abla::qrot(G4double z, G4double a, G4double bet, G4double sig, G4double u, G4double *qr)
00975 {
00976 G4double ucr = 10.0;
00977 G4double dcr = 40.0;
00978 G4double ponq = 0.0, dn = 0.0, n = 0.0, dz = 0.0;
00979
00980 if(((std::fabs(bet)-1.15) < 0) || ((std::fabs(bet)-1.15) == 0)) {
00981 goto qrot10;
00982 }
00983
00984 if((std::fabs(bet)-1.15) > 0) {
00985 goto qrot11;
00986 }
00987
00988 qrot10:
00989 n = a - z;
00990 dz = std::fabs(z - 82.0);
00991 if (n > 104) {
00992 dn = std::fabs(n-126.e0);
00993 }
00994 else {
00995 dn = std::fabs(n - 82.0);
00996 }
00997
00998 bet = 0.022 + 0.003*dn + 0.005*dz;
00999
01000 sig = 25.0*std::pow(bet,2) * sig;
01001
01002 qrot11:
01003 ponq = (u - ucr)/dcr;
01004
01005 if (ponq > 700.0) {
01006 ponq = 700.0;
01007 }
01008 if (sig < 1.0) {
01009 sig = 1.0;
01010 }
01011 (*qr) = 1.0/(1.0 + std::exp(ponq)) * (sig - 1.0) + 1.0;
01012
01013 if ((*qr) < 1.0) {
01014 (*qr) = 1.0;
01015 }
01016
01017 return;
01018 }
01019
01020 void G4Abla::mglw(G4double a, G4double z, G4double *el)
01021 {
01022
01023
01024
01025 G4int a1 = 0, z1 = 0;
01026 G4double xv = 0.0, xs = 0.0, xc = 0.0, xa = 0.0;
01027
01028 a1 = idnint(a);
01029 z1 = idnint(z);
01030
01031 if ((a <= 0.01) || (z < 0.01)) {
01032 (*el) = 1.0e38;
01033 }
01034 else {
01035 xv = -15.56*a;
01036 xs = 17.23*std::pow(a,(2.0/3.0));
01037
01038 if (a > 1.0) {
01039 xc = 0.7*z*(z-1.0)*std::pow((a-1.0),(-1.e0/3.e0));
01040 }
01041 else {
01042 xc = 0.0;
01043 }
01044 }
01045
01046 xa = 23.6*(std::pow((a-2.0*z),2)/a);
01047 (*el) = xv+xs+xc+xa;
01048 return;
01049 }
01050
01051 void G4Abla::mglms(G4double a, G4double z, G4int refopt4, G4double *el)
01052 {
01053
01054
01055
01056
01057
01058
01059
01060
01061
01062
01063
01064
01065
01066
01067
01068
01069
01070
01071 G4int a1 = idnint(a);
01072 G4int z1 = idnint(z);
01073
01074 if ( (a1 <= 0) || (z1 <= 0) || ((a1-z1) <= 0) ) {
01075
01076 (*el) = 0.0;
01077 return;
01078
01079 }
01080 else {
01081
01082
01083 assert(a1 != 0);
01084 (*el) = eflmac(a1,z1,0,refopt4);
01085
01086 if (refopt4 > 0) {
01087 if (refopt4 != 2) {
01088 (*el) = (*el) + ec2sub->ecnz[a1-z1][z1];
01089
01090
01091 }
01092 }
01093 }
01094 return;
01095 }
01096
01097 G4double G4Abla::spdef(G4int a, G4int z, G4int optxfis)
01098 {
01099
01100
01101
01102
01103
01104
01105
01106
01107 G4int index = 0;
01108 G4double x = 0.0, v = 0.0, dx = 0.0;
01109
01110 const G4int alpha2Size = 37;
01111
01112 G4double alpha2[alpha2Size] = {0.0, 2.5464e0, 2.4944e0, 2.4410e0, 2.3915e0, 2.3482e0,
01113 2.3014e0, 2.2479e0, 2.1982e0, 2.1432e0, 2.0807e0, 2.0142e0, 1.9419e0,
01114 1.8714e0, 1.8010e0, 1.7272e0, 1.6473e0, 1.5601e0, 1.4526e0, 1.3164e0,
01115 1.1391e0, 0.9662e0, 0.8295e0, 0.7231e0, 0.6360e0, 0.5615e0, 0.4953e0,
01116 0.4354e0, 0.3799e0, 0.3274e0, 0.2779e0, 0.2298e0, 0.1827e0, 0.1373e0,
01117 0.0901e0, 0.0430e0, 0.0000e0};
01118
01119 dx = 0.02;
01120 x = fissility(a,z,optxfis);
01121
01122 if (x > 1.0) {
01123 x = 1.0;
01124 }
01125
01126 if (x < 0.0) {
01127 x = 0.0;
01128 }
01129
01130 v = (x - 0.3)/dx + 1.0;
01131 index = idnint(v);
01132
01133 if (index < 1) {
01134 return(alpha2[1]);
01135 }
01136
01137 if (index == 36) {
01138 return(alpha2[36]);
01139 }
01140 else {
01141 return(alpha2[index] + (alpha2[index+1] - alpha2[index]) / dx * ( x - (0.3e0 + dx*(index-1))));
01142 }
01143
01144 return alpha2[0];
01145 }
01146
01147 G4double G4Abla::fissility(int a,int z, int optxfis)
01148 {
01149
01150
01151
01152
01153
01154
01155
01156 G4double aa = 0.0, zz = 0.0, i = 0.0;
01157 G4double fissilityResult = 0.0;
01158
01159 aa = double(a);
01160 zz = double(z);
01161 i = double(a-2*z) / aa;
01162
01163
01164 if (optxfis == 0) {
01165 fissilityResult = std::pow(zz,2) / aa /50.8830e0 / (1.0e0 - 1.7826e0 * std::pow(i,2));
01166 }
01167
01168 if (optxfis == 1) {
01169
01170 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));
01171 }
01172
01173 if (optxfis == 2) {
01174
01175 fissilityResult = std::pow(zz,2) / aa /(48.e0*(1.e0 - 17.22e0*std::pow(i,4)));
01176 }
01177
01178 return fissilityResult;
01179 }
01180
01181 void G4Abla::evapora(G4double zprf, G4double aprf, G4double ee, G4double jprf,
01182 G4double *zf_par, G4double *af_par, G4double *mtota_par,
01183 G4double *pleva_par, G4double *pxeva_par, G4double *pyeva_par,
01184 G4int *ff_par, G4int *inttype_par, G4int *inum_par)
01185 {
01186 G4double zf = (*zf_par);
01187 G4double af = (*af_par);
01188 G4double mtota = (*mtota_par);
01189 G4double pleva = (*pleva_par);
01190 G4double pxeva = (*pxeva_par);
01191 G4double pyeva = (*pyeva_par);
01192 G4int ff = (*ff_par);
01193 G4int inttype = (*inttype_par);
01194 G4int inum = (*inum_par);
01195
01196
01197
01198
01199
01200
01201
01202
01203
01204
01205
01206
01207
01208
01209
01210
01211
01212
01213
01214
01215
01216
01217
01218
01219
01220
01221
01222
01223
01224
01225
01226
01227
01228
01229
01230
01231
01232
01233
01234
01235
01236
01237
01238
01239
01240
01241
01242
01243
01244
01245
01246
01247
01248
01249
01250
01251
01252
01253
01254
01255
01256
01257
01258
01259
01260
01261
01262
01263
01264
01265
01266
01267
01268
01269
01270
01271
01272
01273
01274
01275
01276
01277
01278
01279
01280
01281
01282
01283
01284
01285
01286
01287
01288 static G4int sortie = 0;
01289 static G4double epsiln = 0.0, probp = 0.0, probn = 0.0, proba = 0.0, ptotl = 0.0, e = 0.0;
01290 static G4double sn = 0.0, sbp = 0.0, sba = 0.0, x = 0.0, amoins = 0.0, zmoins = 0.0;
01291 G4double ecn = 0.0, ecp = 0.0,eca = 0.0, bp = 0.0, ba = 0.0;
01292 static G4double pteva = 0.0;
01293
01294 static G4int itest = 0;
01295 static G4double probf = 0.0;
01296
01297 static G4int k = 0, j = 0, il = 0;
01298
01299 static G4double ctet1 = 0.0, stet1 = 0.0, phi1 = 0.0;
01300 static G4double sbfis = 0.0, rnd = 0.0;
01301 static G4double selmax = 0.0;
01302 static G4double segs = 0.0;
01303 static G4double ef = 0.0;
01304 static G4int irndm = 0;
01305
01306 static G4double pc = 0.0, malpha = 0.0;
01307
01308 zf = zprf;
01309 af = aprf;
01310 pleva = 0.0;
01311 pteva = 0.0;
01312 pxeva = 0.0;
01313 pyeva = 0.0;
01314
01315 sortie = 0;
01316 ff = 0;
01317
01318 itest = 0;
01319 if (itest == 1) {
01320 G4cout << "***************************" << G4endl;
01321 }
01322
01323 evapora10:
01324
01325 if (itest == 1) {
01326 G4cout <<"------zf,af,ee------" << idnint(zf) << "," << idnint(af) << "," << ee << G4endl;
01327 }
01328
01329
01330
01331 direct(zf,af,ee,jprf,&probp,&probn,&proba,&probf,&ptotl,
01332 &sn,&sbp,&sba,&ecn,&ecp,&eca,&bp,&ba,inttype,inum,itest);
01333
01334
01335
01336
01337 assert((eca+ba) >= 0);
01338 assert((ecp+bp) >= 0);
01339
01340
01341
01342
01343 k = idnint(zf);
01344 j = idnint(af-zf);
01345
01346
01347
01348
01349
01350 il = idnint(jprf);
01351 barfit(k,k+j,il,&sbfis,&segs,&selmax);
01352
01353
01354 if ((fiss->optshp == 1) || (fiss->optshp == 3)) {
01355
01356 fb->efa[j][k] = G4double(sbfis) - ecld->ecgnz[j][k];
01357 }
01358 else {
01359 fb->efa[j][k] = G4double(sbfis);
01360
01361 }
01362 ef = fb->efa[j][k];
01363
01364
01365
01366 if ((sortie == 1) || (ptotl == 0.e0)) {
01367 e = dmin1(sn,sbp,sba);
01368 if (e > 1.0e30) {
01369 if(verboseLevel > 2) {
01370 G4cout <<"erreur a la sortie evapora,e>1.e30,af=" << af <<" zf=" << zf << G4endl;
01371 }
01372 }
01373 if (zf <= 6.0) {
01374 goto evapora100;
01375 }
01376 if (e < 0.0) {
01377 if (sn == e) {
01378 af = af - 1.e0;
01379 }
01380 else if (sbp == e) {
01381 af = af - 1.0;
01382 zf = zf - 1.0;
01383 }
01384 else if (sba == e) {
01385 af = af - 4.0;
01386 zf = zf - 2.0;
01387 }
01388 if (af < 2.5) {
01389 goto evapora100;
01390 }
01391 goto evapora10;
01392 }
01393 goto evapora100;
01394 }
01395 irndm = irndm + 1;
01396
01397
01398
01399
01400
01401 x = double(haz(1))*ptotl;
01402
01403
01404
01405
01406
01407
01408 itest = 0;
01409 if (x < proba) {
01410
01411 if (itest == 1) {
01412 G4cout <<"< alpha evaporation >" << G4endl;
01413 }
01414 amoins = 4.0;
01415 zmoins = 2.0;
01416 epsiln = sba + eca;
01417 assert((std::pow((1.0 + (eca+ba)/3.72834e3),2) - 1.0) >= 0);
01418 pc = std::sqrt(std::pow((1.0 + (eca+ba)/3.72834e3),2) - 1.0) * 3.72834e3;
01419
01420 malpha = 4.0;
01421
01422
01423 volant->iv = volant->iv + 1;
01424 volant->acv[volant->iv] = 4.;
01425 volant->zpcv[volant->iv] = 2.;
01426 volant->pcv[volant->iv] = pc;
01427 }
01428 else if (x < proba+probp) {
01429
01430 if (itest == 1) {
01431 G4cout <<"< proton evaporation >" << G4endl;
01432 }
01433 amoins = 1.0;
01434 zmoins = 1.0;
01435 epsiln = sbp + ecp;
01436 assert((std::pow((1.0 + (ecp + bp)/9.3827e2),2) - 1.0) >= 0);
01437 pc = std::sqrt(std::pow((1.0 + (ecp + bp)/9.3827e2),2) - 1.0) * 9.3827e2;
01438
01439 malpha = 0.0;
01440
01441 volant->iv = volant->iv + 1;
01442 volant->acv[volant->iv] = 1.0;
01443 volant->zpcv[volant->iv] = 1.;
01444 volant->pcv[volant->iv] = pc;
01445 }
01446 else if (x < proba+probp+probn) {
01447
01448 if (itest == 1) {
01449 G4cout <<"< neutron evaporation >" << G4endl;
01450 }
01451 amoins = 1.0;
01452 zmoins = 0.0;
01453 epsiln = sn + ecn;
01454 assert((std::pow((1.0 + (ecn)/9.3956e2),2) - 1.0) >= 0);
01455 pc = std::sqrt(std::pow((1.0 + (ecn)/9.3956e2),2) - 1.0) * 9.3956e2;
01456
01457 malpha = 0.0;
01458
01459
01460 volant->iv = volant->iv + 1;
01461 volant->acv[volant->iv] = 1.;
01462 volant->zpcv[volant->iv] = 0.;
01463 volant->pcv[volant->iv] = pc;
01464 }
01465 else {
01466
01467
01468
01469
01470 if (itest == 1) {
01471 G4cout <<"< fission >" << G4endl;
01472 }
01473 amoins = 0.0;
01474 zmoins = 0.0;
01475 epsiln = ef;
01476
01477 malpha = 0.0;
01478 pc = 0.0;
01479 ff = 1;
01480
01481 }
01482
01483 if (itest == 1) {
01484 G4cout <<"sn,sbp,sba,ef" << sn << "," << sbp << "," << sba <<"," << ef << G4endl;
01485 G4cout <<"probn,probp,proba,probf,ptotl " <<","<< probn <<","<< probp <<","<< proba <<","<< probf <<","<< ptotl << G4endl;
01486 }
01487
01488
01489 af = af - amoins;
01490 zf = zf - zmoins;
01491 ee = ee - epsiln;
01492 if (ee <= 0.01) {
01493 ee = 0.01;
01494 }
01495 mtota = mtota + malpha;
01496
01497 if(ff == 0) {
01498 standardRandom(&rnd,&(hazard->igraine[8]));
01499 ctet1 = 2.0*rnd - 1.0;
01500 standardRandom(&rnd,&(hazard->igraine[4]));
01501 phi1 = rnd*2.0*3.141592654;
01502 stet1 = std::sqrt(1.0 - std::pow(ctet1,2));
01503
01504 volant->xcv[volant->iv] = stet1*std::cos(phi1);
01505 volant->ycv[volant->iv] = stet1*std::sin(phi1);
01506 volant->zcv[volant->iv] = ctet1;
01507 pxeva = pxeva - pc * volant->xcv[volant->iv];
01508 pyeva = pyeva - pc * volant->ycv[volant->iv];
01509 pleva = pleva - pc * ctet1;
01510
01511 }
01512
01513
01514 if ((af < 2.5) || (ff == 1)) {
01515 goto evapora100;
01516 }
01517 goto evapora10;
01518
01519 evapora100:
01520 (*zf_par) = zf;
01521 (*af_par) = af;
01522 (*mtota_par) = mtota;
01523 (*pleva_par) = pleva;
01524 (*pxeva_par) = pxeva;
01525 (*pyeva_par) = pyeva;
01526 (*ff_par) = ff;
01527 (*inttype_par) = inttype;
01528 (*inum_par) = inum;
01529
01530 return;
01531 }
01532
01533 void G4Abla::direct(G4double zprf, G4double a, G4double ee, G4double jprf,
01534 G4double *probp_par, G4double *probn_par, G4double *proba_par,
01535 G4double *probf_par, G4double *ptotl_par, G4double *sn_par,
01536 G4double *sbp_par, G4double *sba_par, G4double *ecn_par,
01537 G4double *ecp_par,G4double *eca_par, G4double *bp_par,
01538 G4double *ba_par, G4int inttype, G4int inum, G4int itest)
01539 {
01540 G4int dummy0 = 0;
01541
01542 G4double probp = (*probp_par);
01543 G4double probn = (*probn_par);
01544 G4double proba = (*proba_par);
01545 G4double probf = (*probf_par);
01546 G4double ptotl = (*ptotl_par);
01547 G4double sn = (*sn_par);
01548 G4double sbp = (*sbp_par);
01549 G4double sba = (*sba_par);
01550 G4double ecn = (*ecn_par);
01551 G4double ecp = (*ecp_par);
01552 G4double eca = (*eca_par);
01553 G4double bp = (*bp_par);
01554 G4double ba = (*ba_par);
01555
01556
01557
01558
01559
01560
01561
01562
01563
01564
01565
01566
01567
01568
01569
01570
01571
01572
01573
01574
01575
01576
01577
01578
01579
01580
01581
01582
01583
01584
01585
01586
01587
01588
01589
01590
01591
01592
01593
01594
01595
01596
01597
01598
01599
01600
01601
01602
01603
01604
01605
01606
01607
01608
01609
01610
01611
01612
01613
01614
01615
01616
01617
01618
01619
01620
01621 static G4double bk = 0.0;
01622 static G4int afp = 0;
01623 static G4double at = 0.0;
01624 static G4double bs = 0.0;
01625 static G4double bshell = 0.0;
01626 static G4double cf = 0.0;
01627 static G4double dconst = 0.0;
01628 static G4double defbet = 0.0;
01629 static G4double denomi = 0.0;
01630 static G4double densa = 0.0;
01631 static G4double densf = 0.0;
01632 static G4double densg = 0.0;
01633 static G4double densn = 0.0;
01634 static G4double densp = 0.0;
01635 static G4double edyn = 0.0;
01636 static G4double eer = 0.0;
01637 static G4double ef = 0.0;
01638 static G4double ft = 0.0;
01639 static G4double ga = 0.0;
01640 static G4double gf = 0.0;
01641 static G4double gn = 0.0;
01642 static G4double gngf = 0.0;
01643 static G4double gp = 0.0;
01644 static G4double gsum = 0.0;
01645 static G4double hbar = 6.582122e-22;
01646 static G4double iflag = 0.0;
01647 static G4int il = 0;
01648 static G4int imaxwell = 0;
01649 static G4int in = 0;
01650 static G4int iz = 0;
01651 static G4int j = 0;
01652 static G4int k = 0;
01653 static G4double ma1z = 0.0;
01654 static G4double ma1z1 = 0.0;
01655 static G4double ma4z2 = 0.0;
01656 static G4double maz = 0.0;
01657 static G4double nprf = 0.0;
01658 static G4double nt = 0.0;
01659 static G4double parc = 0.0;
01660 static G4double pi = 3.14159265;
01661 static G4double pt = 0.0;
01662 static G4double ra = 0.0;
01663 static G4double rat = 0.0;
01664 static G4double refmod = 0.0;
01665 static G4double rf = 0.0;
01666 static G4double rn = 0.0;
01667 static G4double rnd = 0.0;
01668 static G4double rnt = 0.0;
01669 static G4double rp = 0.0;
01670 static G4double rpt = 0.0;
01671 static G4double sa = 0.0;
01672 static G4double sbf = 0.0;
01673 static G4double sbfis = 0.0;
01674 static G4double segs = 0.0;
01675 static G4double selmax = 0.0;
01676 static G4double sp = 0.0;
01677 static G4double tauc = 0.0;
01678 static G4double tconst = 0.0;
01679 static G4double temp = 0.0;
01680 static G4double ts1 = 0.0;
01681 static G4double tsum = 0.0;
01682 static G4double wf = 0.0;
01683 static G4double wfex = 0.0;
01684 static G4double xx = 0.0;
01685 static G4double y = 0.0;
01686
01687 imaxwell = 1;
01688 inttype = 0;
01689
01690
01691
01692 edyn = 1000.0;
01693
01694
01695 if (fiss->bet <= 1.0e-16) {
01696 edyn = 10000.0;
01697 }
01698
01699
01700 eer = ee;
01701 if (inum == 1) {
01702 ilast = 1;
01703 }
01704
01705
01706
01707
01708 refmod = 1;
01709
01710 if (refmod == 1) {
01711 mglms(a,zprf,fiss->optshp,&maz);
01712 mglms(a-1.0,zprf,fiss->optshp,&ma1z);
01713 mglms(a-1.0,zprf-1.0,fiss->optshp,&ma1z1);
01714 mglms(a-4.0,zprf-2.0,fiss->optshp,&ma4z2);
01715 }
01716 else {
01717 mglw(a,zprf,&maz);
01718 mglw(a-1.0,zprf,&ma1z);
01719 mglw(a-1.0,zprf-1.0,&ma1z1);
01720 mglw(a-4.0,zprf-2.0,&ma4z2);
01721 }
01722
01723
01724
01725
01726
01727
01728 sn = ma1z - maz;
01729 sp = ma1z1 - maz;
01730 sa = ma4z2 - maz - 28.29688;
01731 if (zprf < 1.0e0) {
01732 sbp = 1.0e75;
01733 goto direct30;
01734 }
01735
01736
01737
01738
01739 bp = 1.44*(zprf - 1.0)/(2.1*std::pow((a - 1.0),(1.0/3.0)) + 0.0);
01740
01741 sbp = sp + bp;
01742
01743
01744 if (a-4.0 <= 0.0) {
01745 sba = 1.0e+75;
01746 goto direct30;
01747 }
01748
01749
01750
01751
01752 ba = 2.88*(zprf - 2.0)/(2.2*std::pow((a - 4.0),(1.0/3.0)) + 0.0);
01753
01754 sba = sa + ba;
01755
01756
01757 direct30:
01758
01759
01760
01761 if (fiss->ifis > 0) {
01762 k = idnint(zprf);
01763 j = idnint(a - zprf);
01764
01765
01766
01767
01768
01769 il = idnint(jprf);
01770 barfit(k,k+j,il,&sbfis,&segs,&selmax);
01771
01772 if ((fiss->optshp == 1) || (fiss->optshp == 3)) {
01773
01774
01775 fb->efa[j][k] = double(sbfis) - ecld->ecgnz[j][k];
01776 }
01777 else {
01778
01779 fb->efa[j][k] = double(sbfis);
01780 }
01781
01782 ef = fb->efa[j][k];
01783
01784
01785
01786
01787
01788
01789 if (fb->efa[j][k] < 0.0) {
01790 if(verboseLevel > 2) {
01791 G4cout <<"Setting fission barrier to 0" << G4endl;
01792 }
01793 fb->efa[j][k] = 0.0;
01794 }
01795
01796
01797
01798
01799
01800
01801 if (ef < 0.0) {
01802 ef = 0.0;
01803 }
01804 xx = fissility((k+j),k,fiss->optxfis);
01805
01806
01807
01808 y = 1.00 - xx;
01809 if (y < 0.0) {
01810 y = 0.0;
01811 }
01812 if (y > 1.0) {
01813 y = 1.0;
01814 }
01815 bs = bipol(1,y);
01816
01817
01818 bk = bipol(2,y);
01819
01820
01821 }
01822 else {
01823 ef = 1.0e40;
01824 bs = 1.0;
01825 bk = 1.0;
01826 }
01827 sbf = ee - ef;
01828
01829 afp = idnint(a);
01830 iz = idnint(zprf);
01831 in = afp - iz;
01832 bshell = ecld->ecfnz[in][iz];
01833
01834
01835
01836
01837
01838
01839
01840
01841
01842 defbet = 1.58533e0 * spdef(idnint(a),idnint(zprf),fiss->optxfis);
01843
01844
01845
01846
01847
01848
01849
01850
01851
01852
01853
01854 densniv(a,zprf,ee,ef,&densf,bshell,bs,bk,&temp,int(fiss->optshp),int(fiss->optcol),defbet);
01855
01856
01857
01858
01859
01860 ft = temp;
01861 if (iz >= 2) {
01862 bshell = ecld->ecgnz[in][iz-1] - ecld->vgsld[in][iz-1];
01863 defbet = 1.5 * (ecld->alpha[in][iz-1]);
01864
01865
01866 densniv(a-1.0,zprf-1.0e0,ee,sbp,&densp, bshell,1.e0,1.e0,&temp,int(fiss->optshp),int(fiss->optcol),defbet);
01867 assert(temp >= 0);
01868
01869 pt = temp;
01870 if (imaxwell == 1) {
01871
01872
01873
01874 rpt = pt;
01875 ecp = 2.0 * pt;
01876 if(rpt <= 1.0e-3) {
01877 goto direct2914;
01878 }
01879 iflag = 0;
01880 direct1914:
01881 ecp = fmaxhaz(rpt);
01882 iflag = iflag + 1;
01883 if(iflag >= 10) {
01884 standardRandom(&rnd,&(hazard->igraine[5]));
01885 ecp=std::sqrt(rnd)*(eer-sbp);
01886
01887 goto direct2914;
01888 }
01889 if((ecp+sbp) > eer) {
01890 goto direct1914;
01891 }
01892 }
01893 else {
01894 ecp = 2.0 * pt;
01895 }
01896
01897 direct2914:
01898 dummy0 = 0;
01899
01900 }
01901 else {
01902 densp = 0.0;
01903 ecp = 0.0;
01904 pt = 0.0;
01905 }
01906
01907 if (in >= 2) {
01908 bshell = ecld->ecgnz[in-1][iz] - ecld->vgsld[in-1][iz];
01909 defbet = 1.5e0 * (ecld->alpha[in-1][iz]);
01910
01911
01912 densniv(a-1.0,zprf,ee,sn,&densn,bshell, 1.e0,1.e0,&temp,int(fiss->optshp),int(fiss->optcol),defbet);
01913 nt = temp;
01914
01915 if (imaxwell == 1) {
01916
01917
01918
01919 rnt = nt;
01920 ecn = 2.0 * nt;
01921 if(rnt <= 1.e-3) {
01922 goto direct2915;
01923 }
01924
01925 iflag=0;
01926
01927 ecn = fmaxhaz(rnt);
01928 iflag=iflag+1;
01929 if(iflag >= 10) {
01930 standardRandom(&rnd,&(hazard->igraine[6]));
01931 ecn = std::sqrt(rnd)*(eer-sn);
01932
01933 goto direct2915;
01934 }
01935
01936
01937
01938
01939
01940
01941 if((ecn + sn) <= eer) {
01942 ecn = 2.0 * nt;
01943 }
01944 direct2915:
01945 dummy0 = 0;
01946
01947 }
01948 }
01949 else {
01950 densn = 0.0;
01951 ecn = 0.0;
01952 nt = 0.0;
01953 }
01954
01955 if ((in >= 3) && (iz >= 3)) {
01956 bshell = ecld->ecgnz[in-2][iz-2] - ecld->vgsld[in-2][iz-2];
01957 defbet = 1.5 * (ecld->alpha[in-2][iz-2]);
01958
01959
01960 densniv(a-4.0,zprf-2.0e0,ee,sba,&densa,bshell,1.e0,1.e0,&temp,int(fiss->optshp),int(fiss->optcol),defbet);
01961
01962
01963 at = temp;
01964 if (imaxwell == 1) {
01965
01966
01967 rat = at;
01968 eca= 2.e0 * at;
01969 if(rat <= 1.e-3) {
01970 goto direct2916;
01971 }
01972 iflag=0;
01973 direct1916:
01974 eca = fmaxhaz(rat);
01975 iflag=iflag+1;
01976 if(iflag >= 10) {
01977 standardRandom(&rnd,&(hazard->igraine[7]));
01978 eca=std::sqrt(rnd)*(eer-sba);
01979
01980 goto direct2916;
01981 }
01982 if((eca+sba) > eer) {
01983 goto direct1916;
01984 }
01985 else {
01986 eca = 2.0 * at;
01987 }
01988 direct2916:
01989 dummy0 = 0;
01990
01991 }
01992 else {
01993 densa = 0.0;
01994 eca = 0.0;
01995 at = 0.0;
01996 }
01997 }
01998
01999
02000 if (sn < 0.0) {
02001 probn = 1.0;
02002 probp = 0.0;
02003 proba = 0.0;
02004 probf = 0.0;
02005 goto direct70;
02006 }
02007 if (sbp < 0.0) {
02008 probp = 1.0;
02009 probn = 0.0;
02010 proba = 0.0;
02011 probf = 0.0;
02012 goto direct70;
02013 }
02014
02015 if ((a < 50.e0) || (ee > edyn)) {
02016
02017 densf = 0.e0;
02018 }
02019
02020 bshell = ecld->ecgnz[in][iz] - ecld->vgsld[in][iz];
02021 defbet = 1.5e0 * (ecld->alpha[in][iz]);
02022
02023
02024 densniv(a,zprf,ee,0.0e0,&densg,bshell,1.e0,1.e0,&temp,int(fiss->optshp),int(fiss->optcol),defbet);
02025
02026
02027
02028 if ( densg > 0.e0) {
02029
02030
02031 gp = (std::pow(a,(2.0/3.0))/fiss->akap)*densp/densg/pi*std::pow(pt,2);
02032 gn = (std::pow(a,(2.0/3.0))/fiss->akap)*densn/densg/pi*std::pow(nt,2);
02033 ga = (std::pow(a,(2.0/3.0))/fiss->akap)*densa/densg/pi*2.0*std::pow(at,2);
02034 gf = densf/densg/pi/2.0*ft;
02035
02036
02037
02038
02039
02040
02041
02042
02043
02044 if(itest == 1) {
02045 G4cout <<"gn,gp,ga,gf " << gn <<","<< gp <<","<< ga <<","<< gf << G4endl;
02046 }
02047 }
02048 else {
02049 if(verboseLevel > 2) {
02050 G4cout <<"direct: densg <= 0.e0 " << a <<","<< zprf <<","<< ee << G4endl;
02051 }
02052 }
02053
02054 gsum = ga + gp + gn;
02055
02056
02057 if (gsum > 0.0) {
02058 ts1 = hbar / gsum;
02059 }
02060 else {
02061 ts1 = 1.0e99;
02062 }
02063
02064 if (inum > ilast) {
02065 tsum = 0;
02066 }
02067
02068
02069 if (densf == 0.0) {
02070 if (densp == 0.0) {
02071 if (densn == 0.0) {
02072 if (densa == 0.0) {
02073
02074 probf = 0.0;
02075 probp = 0.0;
02076 probn = 0.0;
02077 proba = 0.0;
02078 goto direct70;
02079 }
02080
02081
02082 rf = 0.0;
02083 rp = 0.0;
02084 rn = 0.0;
02085 ra = 1.0;
02086 goto direct50;
02087 }
02088
02089
02090 rf = 0.0;
02091 rp = 0.0;
02092 rn = 1.0;
02093 ra = densa*2.0/densn*std::pow((at/nt),2);
02094 goto direct50;
02095 }
02096
02097 rf = 0.0;
02098 rp = 1.0;
02099 rn = densn/densp*std::pow((nt/pt),2);
02100
02101 ra = densa*2.0/densp*std::pow((at/pt),2);
02102
02103 goto direct50;
02104 }
02105
02106
02107 rf = 1.0;
02108
02109
02110
02111 if (fiss->bet <= 1.0e-16) {
02112 cf = 1.0;
02113 wf = 1.0;
02114 }
02115 else if (sbf > 0.0e0) {
02116 cf = cram(fiss->bet,fiss->homega);
02117
02118
02119
02120 if (ef <= 0.0) {
02121 rp = 0.0;
02122 rn = 0.0;
02123 ra = 0.0;
02124 goto direct50;
02125 }
02126 else {
02127
02128 tauc = tau(fiss->bet,fiss->homega,ef,ft);
02129
02130 }
02131 wfex = (tauc - tsum)/ts1;
02132
02133 if (wfex < 0.0) {
02134 wf = 1.0;
02135 }
02136 else {
02137 wf = std::exp( -wfex);
02138 }
02139 }
02140 else {
02141 cf=1.0;
02142 wf=1.0;
02143 }
02144
02145 if(verboseLevel > 2) {
02146 G4cout <<"tsum,wf,cf " << tsum <<","<< wf <<","<< cf << G4endl;
02147 }
02148
02149 tsum = tsum + ts1;
02150
02151
02152 tconst = 0.7;
02153 dconst = 12.0/std::sqrt(a);
02154
02155 nprf = a - zprf;
02156
02157 if (fiss->optshp >= 2) {
02158 parite(nprf,&parc);
02159
02160 dconst = dconst*parc;
02161 }
02162 else {
02163 dconst= 0.0;
02164 }
02165 if ((ee <= 17.e0) && (fiss->optles == 1) && (iz >= 90) && (in >= 134)) {
02166
02167 gngf = std::pow(a,(2.0/3.0))*tconst/10.0*std::exp((ef-sn+dconst)/tconst);
02168
02169
02170
02171
02172 if (gn == 0.0) {
02173 rn = 0.0;
02174 rp = 0.0;
02175 ra = 0.0;
02176 }
02177 else {
02178 rn=gngf;
02179
02180 rp=gngf*gp/gn;
02181
02182 ra=gngf*ga/gn;
02183
02184 }
02185 } else {
02186
02187
02188
02189
02190 assert(gn > 0 || (gf != 0 && cf != 0));
02191 rn = gn/(gf*cf);
02192
02193
02194
02195 rp = gp/(gf*cf);
02196
02197 ra = ga/(gf*cf);
02198
02199 }
02200 direct50:
02201
02202
02203
02204
02205
02206
02207 denomi = rp+rn+ra+rf;
02208
02209 assert(denomi > 0);
02210
02211 probf = rf/denomi;
02212
02213 probp = rp/denomi;
02214
02215 probn = rn/denomi;
02216
02217 proba = ra/denomi;
02218
02219
02220
02221
02222
02223
02224
02225 assert(std::fabs(probf) <= 1.0);
02226 probf = probf * wf;
02227 if(probf == 1.0) {
02228 probp = 0.0;
02229 probn = 0.0;
02230 proba = 0.0;
02231 }
02232 else {
02233 probp = probp * (wf + (1.e0-wf)/(1.e0-probf));
02234 probn = probn * (wf + (1.e0-wf)/(1.e0-probf));
02235 proba = proba * (wf + (1.e0-wf)/(1.e0-probf));
02236 }
02237 direct70:
02238
02239
02240
02241
02242 ptotl = probp+probn+proba+probf;
02243
02244
02245 ee = eer;
02246 ilast = inum;
02247
02248
02249
02250 (*probp_par) = probp;
02251 (*probn_par) = probn;
02252 (*proba_par) = proba;
02253 (*probf_par) = probf;
02254 (*ptotl_par) = ptotl;
02255 (*sn_par) = sn;
02256 (*sbp_par) = sbp;
02257 (*sba_par) = sba;
02258 (*ecn_par) = ecn;
02259 (*ecp_par) = ecp;
02260 (*eca_par) = eca;
02261 (*bp_par) = bp;
02262 (*ba_par) = ba;
02263 }
02264
02265 void G4Abla::densniv(G4double a, G4double z, G4double ee, G4double esous, G4double *dens, G4double bshell, G4double bs, G4double bk,
02266 G4double *temp, G4int optshp, G4int optcol, G4double defbet)
02267 {
02268
02269
02270
02271
02272
02273
02274
02275
02276
02277
02278
02279
02280
02281
02282
02283
02284
02285
02286
02287
02288
02289
02290
02291
02292
02293
02294
02295
02296
02297
02298
02299
02300
02301
02302
02303
02304
02305
02306
02307
02308
02309
02310
02311
02312
02313
02314
02315
02316 G4double afp = 0.0;
02317 G4double delta0 = 0.0;
02318 G4double deltau = 0.0;
02319 G4double deltpp = 0.0;
02320 G4double e = 0.0;
02321 G4double ecor = 0.0;
02322 G4double ecor1 = 0.0;
02323 G4double ecr = 0.0;
02324 G4double er = 0.0;
02325 G4double fe = 0.0;
02326 G4double fp = 0.0;
02327 G4double he = 0.0;
02328 G4double iz = 0.0;
02329 G4double pa = 0.0;
02330 G4double para = 0.0;
02331 G4double parz = 0.0;
02332 G4double ponfe = 0.0;
02333 G4double ponniv = 0.0;
02334 G4double qr = 0.0;
02335 G4double sig = 0.0;
02336 G4double y01 = 0.0;
02337 G4double y11 = 0.0;
02338 G4double y2 = 0.0;
02339 G4double y21 = 0.0;
02340 G4double y1 = 0.0;
02341 G4double y0 = 0.0;
02342
02343 G4double pi6 = std::pow(3.1415926535,2) / 6.0;
02344 ecr=10.0;
02345 er=28.0;
02346 afp=idnint(a);
02347 iz=idnint(z);
02348
02349
02350 if((ald->optafan == 1)) {
02351 pa = (ald->av)*a + (ald->as)*std::pow(a,(2.e0/3.e0)) + (ald->ak)*std::pow(a,(1.e0/3.e0));
02352 }
02353 else {
02354 pa = (ald->av)*a + (ald->as)*bs*std::pow(a,(2.e0/3.e0)) + (ald->ak)*bk*std::pow(a,(1.e0/3.e0));
02355 }
02356
02357 fp = 0.01377937231e0 * std::pow(a,(5.e0/3.e0)) * (1.e0 + defbet/3.e0);
02358
02359
02360 if (bs > 1.0) {
02361 delta0 = 14;
02362 }
02363 else {
02364 delta0 = 12;
02365 }
02366
02367 if (esous > 1.0e30) {
02368 (*dens) = 0.0;
02369 (*temp) = 0.0;
02370 goto densniv100;
02371 }
02372
02373 e = ee - esous;
02374
02375 if (e < 0.0) {
02376 (*dens) = 0.0;
02377 (*temp) = 0.0;
02378 goto densniv100;
02379 }
02380
02381
02382 if (optshp > 0) {
02383 deltau = bshell;
02384 if (optshp == 2) {
02385 deltau = 0.0;
02386 }
02387 if (optshp >= 2) {
02388
02389
02390 deltpp = -0.25e0* std::pow((delta0/std::sqrt(a)),2) * pa /pi6 + 2.e0*delta0/std::sqrt(a);
02391
02392
02393 parite(a,¶);
02394 if (para < 0.0) {
02395 e = e - delta0/std::sqrt(a);
02396
02397 } else {
02398 parite(z,&parz);
02399 if (parz > 0.e0) {
02400 e = e - 2.0*delta0/std::sqrt(a);
02401
02402 } else {
02403 e = e;
02404
02405 }
02406 }
02407 } else {
02408 deltpp = 0.0;
02409 }
02410 } else {
02411 deltau = 0.0;
02412 deltpp = 0.0;
02413 }
02414 if (e < 0.0) {
02415 e = 0.0;
02416 (*temp) = 0.0;
02417 }
02418
02419
02420 ponfe = -2.5*pa*e*std::pow(a,(-4.0/3.0));
02421
02422 if (ponfe < -700.0) {
02423 ponfe = -700.0;
02424 }
02425 fe = 1.0 - std::exp(ponfe);
02426 if (e < ecr) {
02427
02428 he = 1.0 - std::pow((1.0 - e/ecr),2);
02429 }
02430 else {
02431 he = 1.0;
02432 }
02433
02434
02435
02436 ecor = e + deltau*fe + deltpp*he;
02437
02438 if (ecor <= 0.1) {
02439 ecor = 0.1;
02440 }
02441
02442
02443
02444
02445
02446 if (ee < 5.0) {
02447 y1 = std::sqrt(pa*ecor);
02448
02449 for(int j = 0; j < 5; j++) {
02450 y2 = pa*ecor*(1.e0-std::exp(-y1));
02451
02452 y1 = std::sqrt(y2);
02453
02454 }
02455
02456 y0 = pa/y1;
02457
02458 assert(y0 != 0.0);
02459 (*temp)=1.0/y0;
02460 (*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;
02461 if (ecor < 1.0) {
02462 ecor1=1.0;
02463 y11 = std::sqrt(pa*ecor1);
02464
02465 for(int j = 0; j < 7; j++) {
02466 y21 = pa*ecor1*(1.0-std::exp(-y11));
02467
02468 y11 = std::sqrt(y21);
02469
02470 }
02471
02472 y01 = pa/y11;
02473
02474 (*dens) = (*dens)*std::pow((y01/y0),1.5);
02475 (*temp) = (*temp)*std::pow((y01/y0),1.5);
02476 }
02477 }
02478 else {
02479 ponniv = 2.0*std::sqrt(pa*ecor);
02480
02481 if (ponniv > 700.0) {
02482 ponniv = 700.0;
02483 }
02484
02485
02486 (*dens) = std::pow(pa,(-0.25e0))*std::pow(ecor,(-1.25e0))*std::exp(ponniv) * 0.1477045e0;
02487
02488 (*temp) = std::sqrt(ecor/pa);
02489 }
02490 densniv100:
02491
02492
02493 sig = fp * (*temp);
02494
02495
02496 if (optcol == 1) {
02497 qrot(z,a,defbet,sig,ecor,&qr);
02498 }
02499 else {
02500 qr = 1.0;
02501 }
02502
02503 (*dens) = (*dens) * qr;
02504 }
02505
02506
02507 G4double G4Abla::bfms67(G4double zms, G4double ams)
02508 {
02509
02510
02511
02512
02513
02514 G4double nms = 0.0, ims = 0.0, ksims = 0.0, xms = 0.0, ums = 0.0;
02515
02516 nms = ams - zms;
02517 ims = (nms-zms)/ams;
02518 ksims= 50.15e0 * (1.- 1.78 * std::pow(ims,2));
02519 xms = std::pow(zms,2) / (ams * ksims);
02520 ums = 0.368e0-5.057e0*xms+8.93e0*std::pow(xms,2)-8.71*std::pow(xms,3);
02521 return(0.7322e0*std::pow(zms,2)/std::pow(ams,(0.333333e0))*std::pow(10.e0,ums));
02522 }
02523
02524 void G4Abla::lpoly(G4double x, G4int n, G4double pl[])
02525 {
02526
02527
02528
02529
02530
02531
02532
02533
02534 pl[0] = 1.0;
02535 pl[1] = x;
02536
02537 for(int i = 2; i < n; i++) {
02538 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);
02539 }
02540 }
02541
02542 G4double G4Abla::eflmac(G4int ia, G4int iz, G4int flag, G4int optshp)
02543 {
02544
02545
02546
02547
02548
02549
02550
02551
02552
02553
02554
02555
02556
02557
02558
02559
02560
02561 G4double eflmacResult = 0.0;
02562
02563 G4int in = 0;
02564 G4double z = 0.0, n = 0.0, a = 0.0, av = 0.0, as = 0.0;
02565 G4double a0 = 0.0, c1 = 0.0, c4 = 0.0, b1 = 0.0, b3 = 0.0;
02566 G4double f = 0.0, ca = 0.0, w = 0.0, dp = 0.0, dn = 0.0, dpn = 0.0, efl = 0.0;
02567 G4double rmac = 0.0, bs = 0.0, h = 0.0, r0 = 0.0, kf = 0.0, ks = 0.0;
02568 G4double kv = 0.0, rp = 0.0, ay = 0.0, aden = 0.0, x0 = 0.0, y0 = 0.0;
02569 G4double mh = 0.0, mn = 0.0, esq = 0.0, ael = 0.0, i = 0.0;
02570 G4double pi = 3.141592653589793238e0;
02571
02572
02573
02574 mh = 7.289034;
02575
02576
02577 mn = 8.071431;
02578
02579
02580 esq = 1.4399764;
02581
02582
02583
02584 ael = 1.433e-5;
02585
02586
02587 rp = 0.8;
02588
02589
02590 r0 = 1.16;
02591
02592
02593 ay = 0.68;
02594
02595
02596
02597 aden= 0.70;
02598
02599
02600
02601 rmac= 4.80;
02602
02603
02604 h = 6.6;
02605
02606
02607 w = 30.0;
02608
02609
02610
02611 av = 16.00126;
02612
02613
02614 kv = 1.92240;
02615
02616
02617 as = 21.18466;
02618
02619
02620 ks = 2.345;
02621
02622 a0 = 2.615;
02623
02624
02625 ca = 0.10289;
02626
02627
02628
02629 bs = 1.0;
02630
02631 z = double(iz);
02632 a = double(ia);
02633 in = ia - iz;
02634 n = double(in);
02635 dn = rmac*bs/std::pow(n,(1.0/3.0));
02636 dp = rmac*bs/std::pow(z,(1.0/3.0));
02637 dpn = h/bs/std::pow(a,(2.0/3.0));
02638
02639
02640 c1 = 3.0/5.0*esq/r0;
02641
02642
02643
02644 c4 = 5.0/4.0*std::pow((3.0/(2.0*pi)),(2.0/3.0)) * c1;
02645
02646
02647
02648
02649
02650
02651
02652 kf = std::pow((9.0*pi*z/(4.0*a)),(1.0/3.0))/r0;
02653
02654
02655
02656 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));
02657 i = (n-z)/a;
02658
02659 x0 = r0 * std::pow(a,(1.0/3.0)) / ay;
02660 y0 = r0 * std::pow(a,(1.0/3.0)) / aden;
02661
02662 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);
02663
02664 b3 = 1.0 - 5.0/std::pow(y0,2) * (1.0 - 15.0/(8.0*y0) + 21.0/(8.0 * std::pow(y0,3))
02665 - 3.0/4.0 * (1.0 + 9.0/(2.0*y0) + 7.0/std::pow(y0,2)
02666 + 7.0/(2.0 * std::pow(y0,3))) * std::exp(-2.0*y0));
02667
02668
02669
02670 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
02671 + 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))
02672 + f*std::pow(z,2)/a -ca*(n-z) - ael * std::pow(z,(2.39e0));
02673
02674 if ((in == iz) && (mod(in,2) == 1) && (mod(iz,2) == 1)) {
02675
02676 efl = efl + w*(utilabs(i)+1.e0/a);
02677 }
02678 else {
02679 efl= efl + w* utilabs(i);
02680 }
02681
02682
02683 if (optshp >= 2) {
02684
02685 if ((mod(in,2) == 1) && (mod(iz,2) == 1)) {
02686 efl = efl - dpn;
02687 }
02688 if (mod(in,2) == 1) {
02689 efl = efl + dn;
02690 }
02691 if (mod(iz,2) == 1) {
02692 efl = efl + dp;
02693 }
02694
02695 }
02696
02697 if (flag != 0) {
02698 eflmacResult = (0.5*(dn + dp) - 0.5*dpn);
02699 }
02700 else {
02701 eflmacResult = efl;
02702 }
02703
02704 return eflmacResult;
02705 }
02706
02707 void G4Abla::appariem(G4double a, G4double z, G4double *del)
02708 {
02709
02710
02711
02712
02713
02714 double para = 0.0, parz = 0.0;
02715
02716
02717
02718
02719
02720
02721 parite(a, ¶);
02722
02723 if (para < 0.0) {
02724 (*del) = 0.0;
02725 }
02726 else {
02727 parite(z, &parz);
02728 if (parz > 0.0) {
02729
02730 (*del) = -12.0/std::sqrt(a);
02731 }
02732 else {
02733
02734 (*del) = 12.0/std::sqrt(a);
02735 }
02736 }
02737 }
02738
02739 void G4Abla::parite(G4double n, G4double *par)
02740 {
02741
02742
02743
02744
02745
02746 G4double n1 = 0.0, n2 = 0.0, n3 = 0.0;
02747
02748
02749
02750
02751
02752 n3 = double(idnint(n));
02753 n1 = n3/2.0;
02754 n2 = n1 - dint(n1);
02755
02756 if (n2 > 0.0) {
02757 (*par) = -1.0;
02758 }
02759 else {
02760 (*par) = 1.0;
02761 }
02762 }
02763
02764 G4double G4Abla::tau(G4double bet, G4double homega, G4double ef, G4double t)
02765 {
02766
02767
02768
02769
02770
02771
02772
02773
02774
02775 G4double tauResult = 0.0;
02776
02777 G4double tlim = 8.e0 * ef;
02778 if (t > tlim) {
02779 t = tlim;
02780 }
02781
02782
02783 if (bet/(std::sqrt(2.0)*10.0*(homega/6.582122)) <= 1.0) {
02784 tauResult = std::log(10.0*ef/t)/(bet*1.0e21);
02785
02786 }
02787 else {
02788 tauResult = std::log(10.0*ef/t)/ (2.0*std::pow((10.0*homega/6.582122),2))*(bet*1.0e-21);
02789
02790 }
02791
02792 return tauResult;
02793 }
02794
02795 G4double G4Abla::cram(G4double bet, G4double homega)
02796 {
02797
02798
02799
02800
02801 G4double rel = bet/(20.0*homega/6.582122);
02802 G4double cramResult = std::sqrt(1.0 + std::pow(rel,2)) - rel;
02803
02804
02805 if (cramResult > 1.0) {
02806 cramResult = 1.0;
02807 }
02808
02809
02810 return cramResult;
02811 }
02812
02813 G4double G4Abla::bipol(int iflag, G4double y)
02814 {
02815
02816
02817
02818
02819
02820
02821
02822
02823
02824 int i = 0;
02825
02826 G4double bipolResult = 0.0;
02827
02828 const int bsbkSize = 54;
02829
02830 G4double bk[bsbkSize] = {0.0, 1.00000,1.00087,1.00352,1.00799,1.01433,1.02265,1.03306,
02831 1.04576,1.06099,1.07910,1.10056,1.12603,1.15651,1.19348,
02832 1.23915,1.29590,1.35951,1.41013,1.44103,1.46026,1.47339,
02833 1.48308,1.49068,1.49692,1.50226,1.50694,1.51114,1.51502,
02834 1.51864,1.52208,1.52539,1.52861,1.53177,1.53490,1.53803,
02835 1.54117,1.54473,1.54762,1.55096,1.55440,1.55798,1.56173,
02836 1.56567,1.56980,1.57413,1.57860,1.58301,1.58688,1.58688,
02837 1.58688,1.58740,1.58740, 0.0};
02838
02839 G4double bs[bsbkSize] = {0.0, 1.00000,1.00086,1.00338,1.00750,1.01319,
02840 1.02044,1.02927,1.03974,
02841 1.05195,1.06604,1.08224,1.10085,1.12229,1.14717,1.17623,1.20963,
02842 1.24296,1.26532,1.27619,1.28126,1.28362,1.28458,1.28477,1.28450,
02843 1.28394,1.28320,1.28235,1.28141,1.28042,1.27941,1.27837,1.27732,
02844 1.27627,1.27522,1.27418,1.27314,1.27210,1.27108,1.27006,1.26906,
02845 1.26806,1.26707,1.26610,1.26514,1.26418,1.26325,1.26233,1.26147,
02846 1.26147,1.26147,1.25992,1.25992, 0.0};
02847
02848 i = idint(y/(2.0e-02)) + 1;
02849 assert(i >= 1);
02850
02851 if(i >= bsbkSize) {
02852 if(verboseLevel > 2) {
02853 G4cout <<"G4Abla error: index i = " << i << " is greater than array size permits." << G4endl;
02854 }
02855 bipolResult = 0.0;
02856 }
02857 else {
02858 if (iflag == 1) {
02859 bipolResult = bs[i] + (bs[i+1] - bs[i])/2.0e-02 * (y - 2.0e-02*(i - 1));
02860 }
02861 else {
02862 bipolResult = bk[i] + (bk[i+1] - bk[i])/2.0e-02 * (y - 2.0e-02*(i - 1));
02863 }
02864 }
02865
02866
02867 return bipolResult;
02868 }
02869
02870 void G4Abla::barfit(G4int iz, G4int ia, G4int il, G4double *sbfis, G4double *segs, G4double *selmax)
02871 {
02872
02873
02874
02875
02876
02877
02878
02879
02880
02881
02882
02883
02884
02885
02886
02887
02888
02889
02890
02891
02892
02893
02894
02895
02896
02897
02898
02899
02900
02901
02902
02903
02904
02905
02906
02907
02908
02909
02910
02911
02912
02913
02914
02915
02916
02917
02918
02919
02920
02921
02922
02923
02924
02925
02926
02927
02928
02929
02930
02931
02932
02933
02934
02935
02936
02937
02938 G4double pa[7],pz[7],pl[10];
02939 for(G4int init_i = 0; init_i < 7; init_i++) {
02940 pa[init_i] = 0.0;
02941 pz[init_i] = 0.0;
02942 }
02943 for(G4int init_i = 0; init_i < 10; init_i++) {
02944 pl[init_i] = 0.0;
02945 }
02946
02947 G4double a = 0.0, z = 0.0, amin = 0.0, amax = 0.0, amin2 = 0.0;
02948 G4double amax2 = 0.0, aa = 0.0, zz = 0.0, bfis = 0.0;
02949 G4double bfis0 = 0.0, ell = 0.0, el = 0.0, egs = 0.0, el80 = 0.0, el20 = 0.0;
02950 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;
02951 G4double aj = 0.0, ak = 0.0, a1 = 0.0, a2 = 0.0;
02952
02953 G4int i = 0, j = 0, k = 0, m = 0;
02954 G4int l = 0;
02955
02956 G4double emncof[4][5] = {{-9.01100e+2,-1.40818e+3, 2.77000e+3,-7.06695e+2, 8.89867e+2},
02957 {1.35355e+4,-2.03847e+4, 1.09384e+4,-4.86297e+3,-6.18603e+2},
02958 {-3.26367e+3, 1.62447e+3, 1.36856e+3, 1.31731e+3, 1.53372e+2},
02959 {7.48863e+3,-1.21581e+4, 5.50281e+3,-1.33630e+3, 5.05367e-2}};
02960
02961 G4double elmcof[4][5] = {{1.84542e+3,-5.64002e+3, 5.66730e+3,-3.15150e+3, 9.54160e+2},
02962 {-2.24577e+3, 8.56133e+3,-9.67348e+3, 5.81744e+3,-1.86997e+3},
02963 {2.79772e+3,-8.73073e+3, 9.19706e+3,-4.91900e+3, 1.37283e+3},
02964 {-3.01866e+1, 1.41161e+3,-2.85919e+3, 2.13016e+3,-6.49072e+2}};
02965
02966 G4double emxcof[4][6] = {{9.43596e4,-2.241997e5,2.223237e5,-1.324408e5,4.68922e4,-8.83568e3},
02967 {-1.655827e5,4.062365e5,-4.236128e5,2.66837e5,-9.93242e4,1.90644e4},
02968 {1.705447e5,-4.032e5,3.970312e5,-2.313704e5,7.81147e4,-1.322775e4},
02969 {-9.274555e4,2.278093e5,-2.422225e5,1.55431e5,-5.78742e4,9.97505e3}};
02970
02971 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},
02972 {-1.13269453e+6, 2.97764590e+6,-4.54326326e+6, 3.00464870e+6, -1.44989274e+6,-1.02026610e+5, 6.27959815e+4},
02973 {1.37543304e+6,-3.65808988e+6, 5.47798999e+6,-3.78109283e+6, 1.84131765e+6, 1.53669695e+4,-6.96817834e+4},
02974 {-8.56559835e+5, 2.48872266e+6,-4.07349128e+6, 3.12835899e+6, -1.62394090e+6, 1.19797378e+5, 4.25737058e+4},
02975 {3.28723311e+5,-1.09892175e+6, 2.03997269e+6,-1.77185718e+6, 9.96051545e+5,-1.53305699e+5,-1.12982954e+4},
02976 {4.15850238e+4, 7.29653408e+4,-4.93776346e+5, 6.01254680e+5, -4.01308292e+5, 9.65968391e+4,-3.49596027e+3},
02977 {-1.82751044e+5, 3.91386300e+5,-3.03639248e+5, 1.15782417e+5, -4.24399280e+3,-6.11477247e+3, 3.66982647e+2}};
02978
02979 const G4int sizex = 5;
02980 const G4int sizey = 6;
02981 const G4int sizez = 4;
02982
02983 G4double egscof[sizey][sizey][sizez];
02984
02985 G4double egs1[sizey][sizex] = {{1.927813e5, 7.666859e5, 6.628436e5, 1.586504e5,-7.786476e3},
02986 {-4.499687e5,-1.784644e6,-1.546968e6,-4.020658e5,-3.929522e3},
02987 {4.667741e5, 1.849838e6, 1.641313e6, 5.229787e5, 5.928137e4},
02988 {-3.017927e5,-1.206483e6,-1.124685e6,-4.478641e5,-8.682323e4},
02989 {1.226517e5, 5.015667e5, 5.032605e5, 2.404477e5, 5.603301e4},
02990 {-1.752824e4,-7.411621e4,-7.989019e4,-4.175486e4,-1.024194e4}};
02991
02992 G4double egs2[sizey][sizex] = {{-6.459162e5,-2.903581e6,-3.048551e6,-1.004411e6,-6.558220e4},
02993 {1.469853e6, 6.564615e6, 6.843078e6, 2.280839e6, 1.802023e5},
02994 {-1.435116e6,-6.322470e6,-6.531834e6,-2.298744e6,-2.639612e5},
02995 {8.665296e5, 3.769159e6, 3.899685e6, 1.520520e6, 2.498728e5},
02996 {-3.302885e5,-1.429313e6,-1.512075e6,-6.744828e5,-1.398771e5},
02997 {4.958167e4, 2.178202e5, 2.400617e5, 1.167815e5, 2.663901e4}};
02998
02999 G4double egs3[sizey][sizex] = {{3.117030e5, 1.195474e6, 9.036289e5, 6.876190e4,-6.814556e4},
03000 {-7.394913e5,-2.826468e6,-2.152757e6,-2.459553e5, 1.101414e5},
03001 {7.918994e5, 3.030439e6, 2.412611e6, 5.228065e5, 8.542465e3},
03002 {-5.421004e5,-2.102672e6,-1.813959e6,-6.251700e5,-1.184348e5},
03003 {2.370771e5, 9.459043e5, 9.026235e5, 4.116799e5, 1.001348e5},
03004 {-4.227664e4,-1.738756e5,-1.795906e5,-9.292141e4,-2.397528e4}};
03005
03006 G4double egs4[sizey][sizex] = {{-1.072763e5,-5.973532e5,-6.151814e5, 7.371898e4, 1.255490e5},
03007 {2.298769e5, 1.265001e6, 1.252798e6,-2.306276e5,-2.845824e5},
03008 {-2.093664e5,-1.100874e6,-1.009313e6, 2.705945e5, 2.506562e5},
03009 {1.274613e5, 6.190307e5, 5.262822e5,-1.336039e5,-1.115865e5},
03010 {-5.715764e4,-2.560989e5,-2.228781e5,-3.222789e3, 1.575670e4},
03011 {1.189447e4, 5.161815e4, 4.870290e4, 1.266808e4, 2.069603e3}};
03012
03013 for(i = 0; i < sizey; i++) {
03014 for(j = 0; j < sizex; j++) {
03015
03016
03017
03018
03019 egscof[i][j][0] = egs1[i][j];
03020 egscof[i][j][1] = egs2[i][j];
03021 egscof[i][j][2] = egs3[i][j];
03022 egscof[i][j][3] = egs4[i][j];
03023 }
03024 }
03025
03026
03027 if (iz < 19 || iz > 111) {
03028 goto barfit900;
03029 }
03030
03031 if(iz > 102 && il > 0) {
03032 goto barfit902;
03033 }
03034
03035 z=double(iz);
03036 a=double(ia);
03037 el=double(il);
03038 amin= 1.2e0*z + 0.01e0*z*z;
03039 amax= 5.8e0*z - 0.024e0*z*z;
03040
03041 if(a < amin || a > amax) {
03042 goto barfit910;
03043 }
03044
03045
03046 aa=2.5e-3*a;
03047 zz=1.0e-2*z;
03048 ell=1.0e-2*el;
03049 bfis0 = 0.0;
03050 lpoly(zz,7,pz);
03051 lpoly(aa,7,pa);
03052
03053 for(i = 0; i < 7; i++) {
03054 for(j = 0; j < 7; j++) {
03055 bfis0=bfis0+elzcof[j][i]*pz[i]*pa[j];
03056
03057 }
03058 }
03059
03060 bfis=bfis0;
03061
03062
03063 (*sbfis)=bfis;
03064 egs=0.0;
03065 (*segs)=egs;
03066
03067
03068
03069 amin2 = 1.4e0*z + 0.009e0*z*z;
03070 amax2 = 20.e0 + 3.0e0*z;
03071
03072 if((a < amin2-5.e0 || a > amax2+10.e0) && il > 0) {
03073 goto barfit920;
03074 }
03075
03076 lpoly(zz,5,pz);
03077 lpoly(aa,4,pa);
03078 el80=0.0;
03079 el20=0.0;
03080 elmax=0.0;
03081
03082 for(i = 0; i < 4; i++) {
03083 for(j = 0; j < 5; j++) {
03084
03085
03086 el80 = el80 + elmcof[i][j]*pz[j]*pa[i];
03087 el20 = el20 + emncof[i][j]*pz[j]*pa[i];
03088 }
03089 }
03090
03091 sel80 = el80;
03092 sel20 = el20;
03093
03094
03095 lpoly(zz,6,pz);
03096 lpoly(ell,9,pl);
03097
03098 for(i = 0; i < 4; i++) {
03099 for(j = 0; j < 6; j++) {
03100
03101
03102 elmax = elmax + emxcof[i][j]*pz[j]*pa[i];
03103 }
03104 }
03105
03106
03107 (*selmax)=elmax;
03108
03109
03110 if(il < 1){
03111 return;
03112 }
03113
03114 x = sel20/(*selmax);
03115
03116 y = sel80/(*selmax);
03117
03118
03119 if(el <= sel20) {
03120
03121 q = 0.2/(std::pow(sel20,2)*std::pow(sel80,2)*(sel20-sel80));
03122 qa = q*(4.0*std::pow(sel80,3) - std::pow(sel20,3));
03123 qb = -q*(4.0*std::pow(sel80,2) - std::pow(sel20,2));
03124 bfis = bfis*(1.0 + qa*std::pow(el,2) + qb*std::pow(el,3));
03125 }
03126 else {
03127
03128 aj = (-20.0*std::pow(x,5) + 25.e0*std::pow(x,4) - 4.0)*std::pow((y-1.0),2)*y*y;
03129 ak = (-20.0*std::pow(y,5) + 25.0*std::pow(y,4) - 1.0) * std::pow((x-1.0),2)*x*x;
03130 q = 0.2/(std::pow((y-x)*((1.0-x)*(1.0-y)*x*y),2));
03131 qa = q*(aj*y - ak*x);
03132 qb = -q*(aj*(2.0*y + 1.0) - ak*(2.0*x + 1.0));
03133 z = el/(*selmax);
03134 a1 = 4.0*std::pow(z,5) - 5.0*std::pow(z,4) + 1.0;
03135 a2 = qa*(2.e0*z + 1.e0);
03136 bfis=bfis*(a1 + (z - 1.e0)*(a2 + qb*z)*z*z*(z - 1.e0));
03137 }
03138
03139 if(bfis <= 0.0) {
03140 bfis=0.0;
03141 }
03142
03143 if(el > (*selmax)) {
03144 bfis=0.0;
03145 }
03146 (*sbfis)=bfis;
03147
03148
03149 if(el > (*selmax)) {
03150 return;
03151 }
03152
03153 for(k = 0; k < 4; k++) {
03154 for(l = 0; l < 6; l++) {
03155 for(m = 0; m < 5; m++) {
03156
03157 egs = egs + egscof[l][m][k]*pz[l]*pa[k]*pl[2*m];
03158
03159 }
03160 }
03161 }
03162
03163 (*segs)=egs;
03164 if((*segs) < 0.0) {
03165 (*segs)=0.0;
03166 }
03167
03168 return;
03169
03170 barfit900:
03171 (*sbfis)=0.0;
03172
03173 if (iz < 19) {
03174 (*sbfis) = 1.0e3;
03175 }
03176 (*segs)=0.0;
03177 (*selmax)=0.0;
03178 return;
03179
03180 barfit902:
03181 (*sbfis)=0.0;
03182 (*segs)=0.0;
03183 (*selmax)=0.0;
03184 return;
03185
03186 barfit910:
03187 (*sbfis)=0.0;
03188 (*segs)=0.0;
03189 (*selmax)=0.0;
03190 return;
03191
03192 barfit920:
03193 (*sbfis)=0.0;
03194 (*segs)=0.0;
03195 (*selmax)=0.0;
03196 return;
03197 }
03198
03199 G4double G4Abla::expohaz(G4int k, G4double T)
03200 {
03201
03202
03203
03204 return (-1.0*T*std::log(haz(k)));
03205 }
03206
03207 G4double G4Abla::fd(G4double E)
03208 {
03209
03210
03211 return (E*std::exp(-E));
03212 }
03213
03214 G4double G4Abla::f(G4double E)
03215 {
03216
03217 return (1.0 - (E + 1.0) * std::exp(-E));
03218 }
03219
03220 G4double G4Abla::fmaxhaz(G4double T)
03221 {
03222
03223
03224
03225
03226
03227
03228 const int pSize = 101;
03229 static G4double p[pSize];
03230
03231
03232 static G4int i = 0;
03233 static G4int itest = 0;
03234
03235
03236
03237 p[pSize-1] = 8.0;
03238 G4double x = 0.1;
03239 G4double x1 = 0.0;
03240 G4double y = 0.0;
03241
03242 if (itest == 1) {
03243 goto fmaxhaz120;
03244 }
03245
03246 for(i = 1; i <= 99; i++) {
03247 fmaxhaz20:
03248 x1 = x - (f(x) - double(i)/100.0)/fd(x);
03249 x = x1;
03250 if (std::fabs(f(x) - double(i)/100.0) < 1e-5) {
03251 goto fmaxhaz100;
03252 }
03253 goto fmaxhaz20;
03254 fmaxhaz100:
03255 p[i] = x;
03256 }
03257
03258
03259 itest = 0;
03260
03261
03262 fmaxhaz120:
03263 standardRandom(&y, &(hazard->igraine[17]));
03264 i = nint(y*100);
03265
03266
03267 if(i == 0) {
03268 goto fmaxhaz120;
03269 }
03270
03271 if (i == 1) {
03272 x = p[i]*y*100;
03273 }
03274 else {
03275 x = (p[i] - p[i-1])*(y*100 - i) + p[i];
03276 }
03277
03278 return(x*T);
03279 }
03280
03281 G4double G4Abla::pace2(G4double a, G4double z)
03282 {
03283
03284
03285
03286
03287 G4double pace2 = 0.0;
03288
03289 G4int ii = idint(a+0.5);
03290 G4int jj = idint(z+0.5);
03291
03292 if(ii <= 0 || jj < 0) {
03293 pace2=0.;
03294 return pace2;
03295 }
03296
03297 if(jj > 300) {
03298 pace2=0.0;
03299 }
03300 else {
03301 pace2=pace->dm[ii][jj];
03302 }
03303 pace2=pace2/1000.;
03304
03305 if(pace->dm[ii][jj] == 0.) {
03306 if(ii < 12) {
03307 pace2=-500.;
03308 }
03309 else {
03310 guet(&a, &z, &pace2);
03311 pace2=pace2-ii*931.5;
03312 pace2=pace2/1000.;
03313 }
03314 }
03315
03316 return pace2;
03317 }
03318
03319 void G4Abla::guet(G4double *x_par, G4double *z_par, G4double *find_par)
03320 {
03321
03322
03323
03324
03325
03326
03327 G4double x = (*x_par);
03328 G4double z = (*z_par);
03329 G4double find = (*find_par);
03330
03331 const G4int qrows = 50;
03332 const G4int qcols = 70;
03333 G4double q[qrows][qcols];
03334 for(G4int init_i = 0; init_i < qrows; init_i++) {
03335 for(G4int init_j = 0; init_j < qcols; init_j++) {
03336 q[init_i][init_j] = 0.0;
03337 }
03338 }
03339
03340 G4int ix=G4int(std::floor(x+0.5));
03341 G4int iz=G4int(std::floor(z+0.5));
03342 G4double zz = iz;
03343 G4double xx = ix;
03344 find = 0.0;
03345 G4double avol = 15.776;
03346 G4double asur = -17.22;
03347 G4double ac = -10.24;
03348 G4double azer = 8.0;
03349 G4double xjj = -30.03;
03350 G4double qq = -35.4;
03351 G4double c1 = -0.737;
03352 G4double c2 = 1.28;
03353
03354 if(ix <= 7) {
03355 q[0][1]=939.50;
03356 q[1][1]=938.21;
03357 q[1][2]=1876.1;
03358 q[1][3]=2809.39;
03359 q[2][4]=3728.34;
03360 q[2][3]=2809.4;
03361 q[2][5]=4668.8;
03362 q[2][6]=5606.5;
03363 q[3][5]=4669.1;
03364 q[3][6]=5602.9;
03365 q[3][7]=6535.27;
03366 q[4][6]=5607.3;
03367 q[4][7]=6536.1;
03368 q[5][7]=6548.3;
03369 find=q[iz][ix];
03370 }
03371 else {
03372 G4double xneu=xx-zz;
03373 G4double si=(xneu-zz)/xx;
03374 G4double x13=std::pow(xx,.333);
03375 G4double ee1=c1*zz*zz/x13;
03376 G4double ee2=c2*zz*zz/xx;
03377 G4double aux=1.+(9.*xjj/4./qq/x13);
03378 G4double ee3=xjj*xx*si*si/aux;
03379 G4double ee4=avol*xx+asur*(std::pow(xx,.666))+ac*x13+azer;
03380 G4double tota = ee1 + ee2 + ee3 + ee4;
03381 find = 939.55*xneu+938.77*zz - tota;
03382 }
03383
03384 (*x_par) = x;
03385 (*z_par) = z;
03386 (*find_par) = find;
03387 }
03388
03389
03390
03391
03392 void G4Abla::even_odd(G4double r_origin,G4double r_even_odd,G4int &i_out)
03393 {
03394
03395
03396
03397
03398
03399
03400
03401
03402
03403
03404
03405
03406
03407
03408
03409
03410
03411
03412
03413
03414
03415
03416
03417
03418
03419
03420
03421
03422
03423
03424
03425 G4double r_in = 0.0, r_rest = 0.0, r_help = 0.0;
03426 G4double r_floor = 0.0;
03427 G4double r_middle = 0.0;
03428
03429 G4int n_floor = 0;
03430
03431 r_in = r_origin + 0.5;
03432 r_floor = (float)((int)(r_in));
03433 if (r_even_odd < 0.001) {
03434 i_out = (int)(r_floor);
03435 }
03436 else {
03437 r_rest = r_in - r_floor;
03438 r_middle = r_floor + 0.5;
03439 n_floor = (int)(r_floor);
03440 if (n_floor%2 == 0) {
03441
03442 r_help = r_middle + (r_rest - 0.5) * (1.0 - r_even_odd);
03443 }
03444 else {
03445
03446 r_help = r_middle + (r_rest - 0.5) * (1.0 + r_even_odd);
03447 }
03448 i_out = (int)(r_help);
03449 }
03450 }
03451
03452 G4double G4Abla::umass(G4double z,G4double n,G4double beta)
03453 {
03454
03455
03456
03457
03458
03459
03460
03461
03462 G4double a = 0.0, umass = 0.0;
03463 G4double alpha = 0.0;
03464 G4double xcom = 0.0, xvs = 0.0, xe = 0.0;
03465 const G4double pi = 3.1416;
03466
03467 a = n + z;
03468 alpha = ( std::sqrt(5.0/(4.0*pi)) ) * beta;
03469
03470
03471 xcom = 1.0 - 1.7826 * ((a - 2.0*z)/a)*((a - 2.0*z)/a);
03472
03473
03474 xvs = - xcom * ( 15.4941 * a -
03475 17.9439 * std::pow(a,0.66667) * (1.0+0.4*alpha*alpha) );
03476
03477 xe = z*z * (0.7053/(std::pow(a,0.33333)) * (1.0-0.2*alpha*alpha) - 1.1529/a);
03478
03479 umass = xvs + xe;
03480
03481 return umass;
03482 }
03483
03484 G4double G4Abla::ecoul(G4double z1,G4double n1,G4double beta1,G4double z2,G4double n2,G4double beta2,G4double d)
03485 {
03486
03487
03488
03489
03490
03491
03492
03493
03494
03495
03496
03497
03498
03499
03500 G4double ecoul = 0;
03501 G4double dtot = 0;
03502 const G4double r0 = 1.16;
03503
03504 dtot = r0 * ( std::pow((z1+n1),0.33333) * (1.0+(2.0/3.0)*beta1)
03505 + std::pow((z2+n2),0.33333) * (1.0+(2.0/3.0)*beta2) ) + d;
03506 ecoul = z1 * z2 * 1.44 / dtot;
03507
03508
03509 return ecoul;
03510 }
03511
03512 void G4Abla::fissionDistri(G4double &a,G4double &z,G4double &e,
03513 G4double &a1,G4double &z1,G4double &e1,G4double &v1,
03514 G4double &a2,G4double &z2,G4double &e2,G4double &v2)
03515 {
03516
03517
03518
03519
03520
03521
03522
03523
03524
03525
03526
03527
03528
03529
03530
03531
03532
03533
03534
03535
03536
03537
03538
03539
03540
03541
03542
03543
03544
03545
03546
03547
03548
03549
03550
03551
03552
03553
03554
03555
03556
03557
03558
03559
03560
03561
03562
03563
03564
03565
03566
03567
03568
03569
03570
03571
03572
03573
03574
03575
03576
03577
03578
03579
03580
03581
03582
03583
03584
03585
03586
03587
03588
03589 G4double n = 0.0;
03590 G4double nlight1 = 0.0, nlight2 = 0.0;
03591 G4double aheavy1 = 0.0,alight1 = 0.0, aheavy2 = 0.0, alight2 = 0.0;
03592 G4double eheavy1 = 0.0, elight1 = 0.0, eheavy2 = 0.0, elight2 = 0.0;
03593 G4double zheavy1_shell = 0.0, zheavy2_shell = 0.0;
03594 G4double zlight1 = 0.0, zlight2 = 0.0;
03595 G4double masscurv = 0.0;
03596 G4double sasymm1 = 0.0, sasymm2 = 0.0, ssymm = 0.0, ysum = 0.0, yasymm = 0.0;
03597 G4double ssymm_mode1 = 0.0, ssymm_mode2 = 0.0;
03598 G4double cz_asymm1_saddle = 0.0, cz_asymm2_saddle = 0.0;
03599
03600 G4double wzasymm1_saddle, wzasymm2_saddle, wzsymm_saddle = 0.0;
03601 G4double wzasymm1_scission = 0.0, wzasymm2_scission = 0.0, wzsymm_scission = 0.0;
03602 G4double wzasymm1 = 0.0, wzasymm2 = 0.0, wzsymm = 0.0;
03603 G4double nlight1_eff = 0.0, nlight2_eff = 0.0;
03604 G4int imode = 0;
03605 G4double rmode = 0.0;
03606 G4double z1mean = 0.0, z2mean = 0.0, z1width = 0.0, za1width = 0.0;
03607
03608 G4double n1r = 0.0, n2r = 0.0, a1r = 0.0, a2r = 0.0, n1 = 0.0, n2 = 0.0;
03609
03610 G4double zsymm = 0.0, nsymm = 0.0, asymm = 0.0;
03611 G4double n1mean = 0.0, n2mean, n1width;
03612 G4double dueff = 0.0;
03613
03614 G4double eld = 0.0;
03615
03616 G4double re1 = 0.0, re2 = 0.0, re3 = 0.0;
03617 G4double eps1 = 0.0, eps2 = 0.0;
03618 G4double n1ucd = 0.0, n2ucd = 0.0, z1ucd = 0.0, z2ucd = 0.0;
03619 G4double beta = 0.0, beta1 = 0.0, beta2 = 0.0;
03620
03621 G4double dn1_pol = 0.0;
03622
03623
03624 G4int i_help = 0;
03625
03626
03627 G4double a_levdens = 0.0;
03628
03629 G4double a_levdens_light1 = 0.0, a_levdens_light2 = 0.0;
03630 G4double a_levdens_heavy1 = 0.0, a_levdens_heavy2 = 0.0;
03631 const G4double r_null = 1.16;
03632
03633 G4double epsilon_1_saddle = 0.0, epsilon0_1_saddle = 0.0;
03634 G4double epsilon_2_saddle = 0.0, epsilon0_2_saddle = 0.0, epsilon_symm_saddle = 0.0;
03635 G4double epsilon_1_scission = 0.0, epsilon0_1_scission = 0.0;
03636 G4double epsilon_2_scission = 0.0, epsilon0_2_scission = 0.0;
03637 G4double epsilon_symm_scission = 0.0;
03638
03639 G4double e_eff1_saddle = 0.0, e_eff2_saddle = 0.0;
03640 G4double epot0_mode1_saddle = 0.0, epot0_mode2_saddle = 0.0, epot0_symm_saddle = 0.0;
03641 G4double epot_mode1_saddle = 0.0, epot_mode2_saddle = 0.0, epot_symm_saddle = 0.0;
03642 G4double e_defo = 0.0, e_defo1 = 0.0, e_defo2 = 0.0, e_scission = 0.0, e_asym = 0.0;
03643 G4double e1exc = 0.0, e2exc = 0.0;
03644 G4double e1exc_sigma = 0.0, e2exc_sigma = 0.0;
03645 G4double e1final = 0.0, e2final = 0.0;
03646
03647 const G4double r0 = 1.16;
03648 G4double tker = 0.0;
03649 G4double ekin1 = 0.0, ekin2 = 0.0;
03650
03651 G4double ekinr1 = 0.0, ekinr2 = 0.0;
03652 G4int icz = 0, k = 0;
03653
03654
03655
03656
03657
03658
03659
03660
03661
03662
03663
03664 const G4double nheavy1 = 83.0;
03665
03666 const G4double nheavy2 = 90.0;
03667
03668 const G4double delta_u1_shell = -2.65;
03669
03670
03671 const G4double delta_u2_shell = -3.8;
03672
03673
03674 G4double delta_u1 = 0.0;
03675
03676 G4double delta_u2 = 0.0;
03677
03678 const G4double cz_asymm1_shell = 0.7;
03679
03680 const G4double cz_asymm2_shell = 0.15;
03681
03682 const G4double fwidth_asymm1 = 0.63;
03683
03684 const G4double fwidth_asymm2 = 0.97;
03685
03686
03687 const G4double xlevdens = 12.0;
03688
03689 const G4double fgamma1 = 2.0;
03690
03691 G4double gamma = 0.0;
03692
03693 G4double gamma_heavy1 = 0.0;
03694
03695 G4double gamma_heavy2 = 0.0;
03696
03697 const G4double e_zero_point = 0.5;
03698
03699 G4double e_saddle_scission = 0.0;
03700
03701 const G4double friction_factor = 1.0;
03702
03703
03704
03705 G4double ysymm = 0.0;
03706
03707 G4double yasymm1 = 0.0;
03708
03709 G4double yasymm2 = 0.0;
03710
03711 G4double nheavy1_eff = 0.0;
03712
03713 G4double zheavy1 = 0.0;
03714
03715 G4double nheavy2_eff = 0.0;
03716
03717 G4double zheavy2 = 0.0;
03718
03719 G4double eexc1_saddle = 0.0;
03720
03721 G4double eexc2_saddle = 0.0;
03722
03723 G4double eexc_max = 0.0;
03724
03725 G4double aheavy1_mean = 0.0;
03726
03727 G4double aheavy2_mean = 0.0;
03728
03729 G4double wasymm_saddle = 0.0;
03730
03731 G4double waheavy1_saddle = 0.0;
03732
03733 G4double waheavy2_saddle = 0.0;
03734
03735 G4double wasymm = 0.0;
03736
03737 G4double waheavy1 = 0.0;
03738
03739 G4double waheavy2 = 0.0;
03740
03741 G4double r_e_o = 0.0, r_e_o_exp = 0.0;
03742
03743 G4double cz_symm = 0.0;
03744
03745 G4double cn = 0.0;
03746
03747 G4double cz = 0.0;
03748
03749 const G4double sigzmin = 1.16;
03750
03751 const G4double d = 2.0;
03752
03753
03754
03755 const G4double cpol1 = 0.65;
03756
03757 const G4double cpol2 = 0.55;
03758
03759 const G4int nzpol = 1;
03760
03761 const G4int itest = 0;
03762
03763
03764 G4double reps1 = 0.0, reps2 = 0.0, rn1_pol = 0.0;
03765
03766 G4int kkk = 0;
03767
03768
03769
03770
03771 if(itest == 1) {
03772 G4cout << " cn mass " << a << G4endl;
03773 G4cout << " cn charge " << z << G4endl;
03774 G4cout << " cn energy " << e << G4endl;
03775 }
03776
03777
03778 n = a - z;
03779
03780 k = 0;
03781 icz = 0;
03782 if ( (std::pow(z,2)/a < 25.0) || (n < nheavy2) || (e > 500.0) ) {
03783 icz = -1;
03784
03785 goto milledeux;
03786 }
03787
03788 nlight1 = n - nheavy1;
03789 nlight2 = n - nheavy2;
03790
03791
03792
03793
03794
03795
03796 zheavy1_shell = ((nheavy1/n) * z) - ((a/n) * cpol1);
03797 zheavy2_shell = ((nheavy2/n) * z) - ((a/n) * cpol2);
03798
03799 e_saddle_scission =
03800 (-24.0 + 0.02227 * (std::pow(z,2))/(std::pow(a,0.33333)) ) * friction_factor;
03801
03802
03803
03804
03805 if (e_saddle_scission > 0.) {
03806 e_saddle_scission = e_saddle_scission;
03807 }
03808 else {
03809 e_saddle_scission = 0.;
03810 }
03811
03812
03813
03814
03815
03816
03817
03818
03819
03820
03821
03822 if ( (std::pow(z,2))/a < 34.0) {
03823 masscurv = std::pow( 10.0,(-1.093364 + 0.082933 * (std::pow(z,2)/a)
03824 - 0.0002602 * (std::pow(z,4)/std::pow(a,2))) );
03825 } else {
03826 masscurv = std::pow( 10.0,(3.053536 - 0.056477 * (std::pow(z,2)/a)
03827 + 0.0002454 * (std::pow(z,4)/std::pow(a,2))) );
03828 }
03829
03830 cz_symm = (8.0/std::pow(z,2)) * masscurv;
03831
03832 if(itest == 1) {
03833 G4cout << "cz_symmetry= " << cz_symm << G4endl;
03834 }
03835
03836 if (cz_symm < 0) {
03837 icz = -1;
03838
03839 goto milledeux;
03840 }
03841
03842
03843 zsymm = z/2.0;
03844 nsymm = n/2.0;
03845 asymm = nsymm + zsymm;
03846
03847 zheavy1 = (cz_symm*zsymm + cz_asymm1_shell*zheavy1_shell)/(cz_symm + cz_asymm1_shell);
03848 zheavy2 = (cz_symm*zsymm + cz_asymm2_shell*zheavy2_shell)/(cz_symm + cz_asymm2_shell);
03849
03850 nheavy1_eff = (zheavy1 + (a/n * cpol1))*(n/z);
03851 nheavy2_eff = (zheavy2 + (a/n * cpol2))*(n/z);
03852 nlight1_eff = n - nheavy1_eff;
03853 nlight2_eff = n - nheavy2_eff;
03854
03855 zlight1 = z - zheavy1;
03856
03857 zlight2 = z - zheavy2;
03858 aheavy1 = nheavy1_eff + zheavy1;
03859 aheavy2 = nheavy2_eff + zheavy2;
03860 aheavy1_mean = aheavy1;
03861 aheavy2_mean = aheavy2;
03862 alight1 = nlight1_eff + zlight1;
03863 alight2 = nlight2_eff + zlight2;
03864
03865 a_levdens = a / xlevdens;
03866 a_levdens_heavy1 = aheavy1 / xlevdens;
03867 a_levdens_heavy2 = aheavy2 / xlevdens;
03868 a_levdens_light1 = alight1 / xlevdens;
03869 a_levdens_light2 = alight2 / xlevdens;
03870 gamma = a_levdens / (0.4 * (std::pow(a,1.3333)) );
03871 gamma_heavy1 = ( a_levdens_heavy1 / (0.4 * (std::pow(aheavy1,1.3333)) ) ) * fgamma1;
03872 gamma_heavy2 = a_levdens_heavy2 / (0.4 * (std::pow(aheavy2,1.3333)) );
03873
03874 cz_asymm1_saddle = cz_asymm1_shell + cz_symm;
03875 cz_asymm2_saddle = cz_asymm2_shell + cz_symm;
03876
03877
03878
03879 cn = umass(zsymm,(nsymm+1.),0.0) + umass(zsymm,(nsymm-1.),0.0)
03880 + 1.44 * (std::pow(zsymm,2))/
03881 ( (std::pow(r_null,2)) *
03882 ( std::pow((asymm+1.0),0.33333) + std::pow((asymm-1.0),0.33333) ) *
03883 ( std::pow((asymm+1.0),0.33333) + std::pow((asymm-1.0),0.33333) ) )
03884 - 2.0 * umass(zsymm,nsymm,0.0)
03885 - 1.44 * (std::pow(zsymm,2))/
03886 ( ( 2.0 * r_null * (std::pow(asymm,0.33333)) ) *
03887 ( 2.0 * r_null * (std::pow(asymm,0.33333)) ) );
03888
03889
03890 delta_u1 = delta_u1_shell + (std::pow((zheavy1_shell-zheavy1),2))*cz_asymm1_shell;
03891
03892 delta_u2 = delta_u2_shell + (std::pow((zheavy2_shell-zheavy2),2))*cz_asymm2_shell;
03893
03894
03895
03896
03897 epot0_mode1_saddle = (std::pow((zheavy1-zsymm),2)) * cz_symm;
03898 epot0_mode2_saddle = (std::pow((zheavy2-zsymm),2)) * cz_symm;
03899 epot0_symm_saddle = 0.0;
03900
03901 if (itest == 1) {
03902 G4cout << "check zheavy1 = " << zheavy1 << G4endl;
03903 G4cout << "check zheavy2 = " << zheavy2 << G4endl;
03904 G4cout << "check zsymm = " << zsymm << G4endl;
03905 G4cout << "check czsymm = " << cz_symm << G4endl;
03906 G4cout << "check epot0_mode1_saddle = " << epot0_mode1_saddle << G4endl;
03907 G4cout << "check epot0_mode2_saddle = " << epot0_mode2_saddle << G4endl;
03908 G4cout << "check epot0_symm_saddle = " << epot0_symm_saddle << G4endl;
03909 G4cout << "delta_u1 = " << delta_u1 << G4endl;
03910 G4cout << "delta_u2 = " << delta_u2 << G4endl;
03911 }
03912
03913
03914
03915
03916 epot_mode1_saddle = epot0_mode1_saddle + delta_u1;
03917 epot_mode2_saddle = epot0_mode2_saddle + delta_u2;
03918 epot_symm_saddle = epot0_symm_saddle;
03919 if (itest == 1) {
03920 G4cout << "check epot_mode1_saddle = " << epot_mode1_saddle << G4endl;
03921 G4cout << "check epot_mode2_saddle = " << epot_mode2_saddle << G4endl;
03922 G4cout << "check epot_symm_saddle = " << epot_symm_saddle << G4endl;
03923 }
03924
03925
03926 dueff = min(epot_mode1_saddle,epot_mode2_saddle);
03927 dueff = min(dueff,epot_symm_saddle);
03928 dueff = dueff - epot_symm_saddle;
03929
03930 eld = e + dueff + e_zero_point;
03931
03932 if (itest == 1) {
03933 G4cout << "check dueff = " << dueff << G4endl;
03934 G4cout << "check e = " << e << G4endl;
03935 G4cout << "check e_zero_point = " << e_zero_point << G4endl;
03936 G4cout << "check eld = " << eld << G4endl;
03937 }
03938
03939
03940
03941
03942
03943
03944
03945
03946
03947 eheavy1 = e * aheavy1 / a;
03948 eheavy2 = e * aheavy2 / a;
03949 elight1 = e * alight1 / a;
03950 elight2 = e * alight2 / a;
03951
03952 epsilon0_1_saddle = eld - e_zero_point - epot0_mode1_saddle;
03953
03954 epsilon0_2_saddle = eld - e_zero_point - epot0_mode2_saddle;
03955
03956
03957 epsilon_1_saddle = eld - e_zero_point - epot_mode1_saddle;
03958
03959 epsilon_2_saddle = eld - e_zero_point - epot_mode2_saddle;
03960
03961 epsilon_symm_saddle = eld - e_zero_point - epot_symm_saddle;
03962
03963
03964 eexc1_saddle = epsilon_1_saddle;
03965 eexc2_saddle = epsilon_2_saddle;
03966 eexc_max = max(eexc1_saddle,eexc2_saddle);
03967 eexc_max = max(eexc_max,eld);
03968
03969
03970
03971
03972 epsilon0_1_scission = eld + e_saddle_scission - epot0_mode1_saddle;
03973
03974 epsilon0_2_scission = eld + e_saddle_scission - epot0_mode2_saddle;
03975
03976
03977 epsilon_1_scission = eld + e_saddle_scission - epot_mode1_saddle;
03978
03979 epsilon_2_scission = eld+ e_saddle_scission - epot_mode2_saddle;
03980
03981 epsilon_symm_scission = eld + e_saddle_scission - epot_symm_saddle;
03982
03983
03984
03985
03986 e_eff1_saddle = epsilon0_1_saddle - delta_u1 * (std::exp((-epsilon_1_saddle*gamma)));
03987
03988 if (e_eff1_saddle > 0.0) {
03989 wzasymm1_saddle = std::sqrt( (0.5 *
03990 (std::sqrt(1.0/a_levdens*e_eff1_saddle)) /
03991 (cz_asymm1_shell * std::exp((-epsilon_1_saddle*gamma)) + cz_symm) ) );
03992 }
03993 else {
03994 wzasymm1_saddle = 1.0;
03995 }
03996
03997 e_eff2_saddle = epsilon0_2_saddle - delta_u2 * (std::exp((-epsilon_2_saddle*gamma)));
03998 if (e_eff2_saddle > 0.0) {
03999 wzasymm2_saddle = std::sqrt( (0.5 *
04000 (std::sqrt(1.0/a_levdens*e_eff2_saddle)) /
04001 (cz_asymm2_shell * std::exp((-epsilon_2_saddle*gamma)) + cz_symm) ) );
04002 }
04003 else {
04004 wzasymm2_saddle = 1.0;
04005 }
04006
04007 if (eld > e_zero_point) {
04008 if ( (eld + epsilon_symm_saddle) < 0.0) {
04009 G4cout << "<e> eld + epsilon_symm_saddle < 0" << G4endl;
04010 }
04011 wzsymm_saddle = std::sqrt( (0.5 *
04012 (std::sqrt(1.0/a_levdens*(eld+epsilon_symm_saddle))) / cz_symm ) );
04013 } else {
04014 wzsymm_saddle = 1.0;
04015 }
04016
04017 if (itest == 1) {
04018 G4cout << "wz1(saddle) = " << wzasymm1_saddle << G4endl;
04019 G4cout << "wz2(saddle) = " << wzasymm2_saddle << G4endl;
04020 G4cout << "wzsymm(saddle) = " << wzsymm_saddle << G4endl;
04021 }
04022
04023
04024
04025
04026 wzsymm_scission = wzsymm_saddle;
04027
04028 if (e_saddle_scission == 0.0) {
04029
04030 wzasymm1_scission = wzasymm1_saddle;
04031 wzasymm2_scission = wzasymm2_saddle;
04032
04033 }
04034 else {
04035
04036 if (nheavy1_eff > 75.0) {
04037 wzasymm1_scission = (std::sqrt(21.0)) * z/a;
04038 wzasymm2_scission = (std::sqrt (max( (70.0-28.0)/3.0*(z*z/a-35.0)+28.,0.0 )) ) * z/a;
04039 }
04040 else {
04041 wzasymm1_scission = wzasymm1_saddle;
04042 wzasymm2_scission = wzasymm2_saddle;
04043 }
04044
04045 }
04046
04047 wzasymm1_scission = max(wzasymm1_scission,wzasymm1_saddle);
04048 wzasymm2_scission = max(wzasymm2_scission,wzasymm2_saddle);
04049
04050 wzasymm1 = wzasymm1_scission * fwidth_asymm1;
04051 wzasymm2 = wzasymm2_scission * fwidth_asymm2;
04052 wzsymm = wzsymm_scission;
04053
04054
04055
04056
04057
04058
04059
04060
04061
04062
04063
04064
04065 wasymm = wzsymm * a/z;
04066 waheavy1 = wzasymm1 * a/z;
04067 waheavy2 = wzasymm2 * a/z;
04068
04069 wasymm_saddle = wzsymm_saddle * a/z;
04070 waheavy1_saddle = wzasymm1_saddle * a/z;
04071 waheavy2_saddle = wzasymm2_saddle * a/z;
04072
04073 if (itest == 1) {
04074 G4cout << "wasymm = " << wzsymm << G4endl;
04075 G4cout << "waheavy1 = " << waheavy1 << G4endl;
04076 G4cout << "waheavy2 = " << waheavy2 << G4endl;
04077 }
04078
04079
04080 if ( (epsilon0_1_saddle - delta_u1*std::exp((-epsilon_1_saddle*gamma_heavy1))) < 0.0) {
04081 sasymm1 = -10.0;
04082 }
04083 else {
04084 sasymm1 = 2.0 * std::sqrt( a_levdens * (epsilon0_1_saddle -
04085 delta_u1*(std::exp((-epsilon_1_saddle*gamma_heavy1))) ) );
04086 }
04087
04088 if ( (epsilon0_2_saddle - delta_u2*std::exp((-epsilon_2_saddle*gamma_heavy2))) < 0.0) {
04089 sasymm2 = -10.0;
04090 }
04091 else {
04092 sasymm2 = 2.0 * std::sqrt( a_levdens * (epsilon0_2_saddle -
04093 delta_u2*(std::exp((-epsilon_2_saddle*gamma_heavy2))) ) );
04094 }
04095
04096 if (epsilon_symm_saddle > 0.0) {
04097 ssymm = 2.0 * std::sqrt( a_levdens*(epsilon_symm_saddle) );
04098 }
04099 else {
04100 ssymm = -10.0;
04101 }
04102
04103 if (ssymm > -10.0) {
04104 ysymm = 1.0;
04105
04106 if (epsilon0_1_saddle < 0.0) {
04107
04108 yasymm1 = std::exp((sasymm1-ssymm)) * wzasymm1_saddle / wzsymm_saddle * 2.0;
04109
04110 }
04111 else {
04112
04113 ssymm_mode1 = 2.0 * std::sqrt( a_levdens*(epsilon0_1_saddle) );
04114 yasymm1 = ( std::exp((sasymm1-ssymm)) - std::exp((ssymm_mode1 - ssymm)) )
04115 * wzasymm1_saddle / wzsymm_saddle * 2.0;
04116 }
04117
04118 if (epsilon0_2_saddle < 0.0) {
04119
04120 yasymm2 = std::exp((sasymm2-ssymm)) * wzasymm2_saddle / wzsymm_saddle * 2.0;
04121
04122 }
04123 else {
04124
04125 ssymm_mode2 = 2.0 * std::sqrt( a_levdens*(epsilon0_2_saddle) );
04126 yasymm2 = ( std::exp((sasymm2-ssymm)) - std::exp((ssymm_mode2 - ssymm)) )
04127 * wzasymm2_saddle / wzsymm_saddle * 2.0;
04128 }
04129
04130
04131
04132 }
04133 else {
04134 if ( (sasymm1 > -10.0) && (sasymm2 > -10.0) ) {
04135 ysymm = 0.0;
04136 yasymm1 = std::exp(sasymm1) * wzasymm1_saddle * 2.0;
04137 yasymm2 = std::exp(sasymm2) * wzasymm2_saddle * 2.0;
04138 }
04139 }
04140
04141
04142 ysum = ysymm + yasymm1 + yasymm2;
04143 if (ysum > 0.0) {
04144 ysymm = ysymm / ysum;
04145 yasymm1 = yasymm1 / ysum;
04146 yasymm2 = yasymm2 / ysum;
04147 yasymm = yasymm1 + yasymm2;
04148 }
04149 else {
04150 ysymm = 0.0;
04151 yasymm1 = 0.0;
04152 yasymm2 = 0.0;
04153
04154 if ( (epsilon_symm_saddle < epsilon_1_saddle) && (epsilon_symm_saddle < epsilon_2_saddle) ) {
04155 ysymm = 1.0;
04156 }
04157 else {
04158 if (epsilon_1_saddle < epsilon_2_saddle) {
04159 yasymm1 = 1.0;
04160 }
04161 else {
04162 yasymm2 = 1.0;
04163 }
04164 }
04165 }
04166
04167 if (itest == 1) {
04168 G4cout << "ysymm normalized= " << ysymm << G4endl;
04169 G4cout << "yasymm1 normalized= " << yasymm1 << G4endl;
04170 G4cout << "yasymm2 normalized= " << yasymm2 << G4endl;
04171 }
04172
04173
04174
04175
04176 if ((int)(z) % 2 == 0) {
04177 r_e_o_exp = -0.017 * (e_saddle_scission + eld) * (e_saddle_scission + eld);
04178 if ( r_e_o_exp < -307.0) {
04179 r_e_o_exp = -307.0;
04180 r_e_o = std::pow(10.0,r_e_o_exp);
04181 }
04182 else {
04183 r_e_o = std::pow(10.0,r_e_o_exp);
04184 }
04185 }
04186 else {
04187 r_e_o = 0.0;
04188 }
04189
04190
04191
04192
04193
04194
04195
04196
04197
04198
04199
04200
04201
04202 do {
04203 rmode = haz(k);
04204
04205
04206
04207 if (rmode < yasymm1) {
04208 imode = 1;
04209 } else if ( (rmode > yasymm1) && (rmode < (yasymm1+yasymm2)) ) {
04210 imode = 2;
04211 } else if ( (rmode > yasymm1) && (rmode > (yasymm1+yasymm2)) ) {
04212 imode = 3;
04213 }
04214 } while(imode == 0);
04215
04216
04217
04218
04219 if (imode == 1) {
04220 z1mean = zheavy1;
04221 z1width = wzasymm1;
04222 }
04223 if (imode == 2) {
04224 z1mean = zheavy2;
04225 z1width = wzasymm2;
04226 }
04227 if (imode == 3) {
04228 z1mean = zsymm;
04229 z1width = wzsymm;
04230 }
04231
04232 if (itest == 1) {
04233 G4cout << "nbre aleatoire tire " << rmode << G4endl;
04234 G4cout << "fission mode " << imode << G4endl;
04235 G4cout << "z1mean= " << z1mean << G4endl;
04236 G4cout << "z1width= " << z1width << G4endl;
04237 }
04238
04239
04240 z1 = 1.0;
04241 z2 = 1.0;
04242 while ( (z1<5.0) || (z2<5.0) ) {
04243
04244
04245 z1 = gausshaz(k, z1mean, z1width);
04246 z2 = z - z1;
04247 }
04248 if (itest == 1) {
04249 G4cout << "ff charge sample " << G4endl;
04250 G4cout << "z1 = " << z1 << G4endl;
04251 G4cout << "z2 = " << z2 << G4endl;
04252 }
04253
04254
04255
04256
04257
04258 z2 = z - z1;
04259
04260
04261 if (imode == 1) {
04262 n1mean = (z1 + cpol1 * a/n) * n/z;
04263 }
04264 if (imode == 2) {
04265 n1mean = (z1 + cpol2 * a/n) * n/z;
04266 }
04267
04268
04269
04270
04271
04272
04273
04274
04275
04276
04277
04278
04279
04280
04281
04282
04283 if (imode == 3) {
04284 n1ucd = z1 * n/z;
04285 n2ucd = z2 * n/z;
04286 re1 = umass(z1,n1ucd,0.6) + umass(z2,n2ucd,0.6) + ecoul(z1,n1ucd,0.6,z2,n2ucd,0.6,d);
04287 re2 = umass(z1,n1ucd+1.,0.6) + umass(z2,n2ucd-1.,0.6) + ecoul(z1,n1ucd+1.,0.6,z2,n2ucd-1.,0.6,d);
04288 re3 = umass(z1,n1ucd+2.,0.6) + umass(z2,n2ucd-2.,0.6) + ecoul(z1,n1ucd+2.,0.6,z2,n2ucd-2.,0.6,d);
04289 eps2 = (re1-2.0*re2+re3) / 2.0;
04290 eps1 = re2 - re1 - eps2;
04291 dn1_pol = - eps1 / (2.0 * eps2);
04292 n1mean = n1ucd + dn1_pol;
04293 }
04294
04295 n2mean = n - n1mean;
04296 z2mean = z - z1mean;
04297
04298
04299
04300
04301
04302
04303
04304
04305
04306
04307
04308 n1r = n1mean;
04309 n2r = n2mean;
04310 a1r = n1r + z1;
04311 a2r = n2r + z2;
04312
04313 if (imode == 1) { ;
04315 e_scission = max(epsilon_1_scission,1.0);
04316 if (n1mean > (n * 0.5) ) {
04318 beta1 = 0.0;
04319 beta2 = 0.6;
04320 e1exc = epsilon_1_scission * a1r / a;
04321 e_defo = umass(z2,n2r,beta2) - umass(z2,n2r,0.0);
04322 e2exc = epsilon_1_scission * a2r / a + e_defo;
04323 }
04324 else {
04326 beta1 = 0.6;
04327 beta2 = 0.0;
04328 e_defo = umass(z1,n1r,beta1) - umass(z1,n1r,0.0);
04329 e1exc = epsilon_1_scission * a1r / a + e_defo;
04330 e2exc = epsilon_1_scission * a2r / a;
04331 }
04332 }
04333
04334 if (imode == 2) {
04336 e_scission = max(epsilon_2_scission,1.0);
04337 if (n1mean > (n * 0.5) ) {
04339 beta1 = (n1r - nheavy2) * 0.034 + 0.3;
04340 e_defo = umass(z1,n1r,beta1) - umass(z1,n1r,0.0);
04341 e1exc = epsilon_2_scission * a1r / a + e_defo;
04342 beta2 = 0.6 - beta1;
04343 e_defo = umass(z2,n2r,beta2) - umass(z2,n2r,0.0);
04344 e2exc = epsilon_2_scission * a2r / a + e_defo;
04345 }
04346 else {
04348 beta2 = (n2r - nheavy2) * 0.034 + 0.3;
04349 e_defo = umass(z2,n2r,beta2) - umass(z2,n2r,0.0);
04350 e2exc = epsilon_2_scission * a2r / a + e_defo;
04351 beta1 = 0.6 - beta2;
04352 e_defo = umass(z1,n1r,beta1) - umass(z1,n1r,0.0);
04353 e1exc = epsilon_2_scission * a1r / a + e_defo;
04354 }
04355 }
04356
04357 if (imode == 3) {
04358
04359
04360
04361
04362
04363
04364
04365 beta = 0.177963 + 0.0153241 * zsymm - 0.000162037 * zsymm*zsymm;
04366 beta1 = 0.177963 + 0.0153241 * z1 - 0.000162037 * z1*z1;
04367
04368 e_defo1 = umass(z1,n1r,beta1) - umass(z1,n1r,0.0);
04369 beta2 = 0.177963 + 0.0153241 * z2 - 0.000162037 * z2*z2;
04370
04371 e_defo2 = umass(z2,n2r,beta2) - umass(z2,n2r,0.0);
04372 e_asym = umass(z1 , n1r, beta1) + umass(z2, n2r ,beta2)
04373 + ecoul(z1,n1r,beta1,z2,n2r,beta2,2.0)
04374 - 2.0 * umass(zsymm,nsymm,beta)
04375 - ecoul(zsymm,nsymm,beta,zsymm,nsymm,beta,2.0);
04376
04377 e_scission = max((epsilon_symm_scission - e_asym),1.0);
04378
04379 e1exc = e_scission * a1r / a + e_defo1;
04380 e2exc = e_scission * a2r / a + e_defo2;
04381 }
04382
04383
04384
04385
04386
04387
04388
04389
04390
04391
04392
04393
04394
04395
04396
04397
04398
04399
04400 if ( (imode == 1) || (imode == 2) ) {
04401 cn=(umass(z1,n1mean+1.,beta1) + umass(z1,n1mean-1.,beta1)
04402 + umass(z2,n2mean+1.,beta2) + umass(z2,n2mean-1.,beta2)
04403 + ecoul(z1,n1mean+1.,beta1,z2,n2mean-1.,beta2,2.0)
04404 + ecoul(z1,n1mean-1.,beta1,z2,n2mean+1.,beta2,2.0)
04405 - 2.0 * ecoul(z1,n1mean,beta1,z2,n2mean,beta2,2.0)
04406 - 2.0 * umass(z1, n1mean, beta1)
04407 - 2.0 * umass(z2, n2mean, beta2) ) * 0.5;
04408
04409
04410
04411
04412
04413
04414
04415 n1width=std::sqrt( (0.5 * (std::sqrt(1.0/a_levdens*(e_scission)))/cn) );
04416 n1width=max(n1width, sigzmin);
04417
04418
04419 n1r = 1.0;
04420 n2r = 1.0;
04421 while ( (n1r<5.0) || (n2r<5.0) ) {
04422
04423
04424 n1r = gausshaz(k, n1mean, n1width);
04425 n2r = n - n1r;
04426 }
04427
04428 if (itest == 1) {
04429 G4cout << "after neutron sample " << n1r << G4endl;
04430 }
04431 n1r = (float)( (int)((n1r+0.5)) );
04432 n2r = n - n1r;
04433
04434 even_odd(z1,r_e_o,i_help);
04435
04436 z1 = (float)(i_help);
04437 z2 = z - z1;
04438
04439 a1r = z1 + n1r;
04440 a2r = z2 + n2r;
04441 }
04442
04443 if (imode == 3) {
04445 if (nzpol > 0.0) {
04446
04447 cz = ( umass(z1-1., n1mean+1.,beta1)
04448 + umass(z2+1., n2mean-1.,beta1)
04449 + umass(z1+1., n1mean-1.,beta2)
04450 + umass(z2 - 1., n2mean + 1.,beta2)
04451 + ecoul(z1-1.,n1mean+1.,beta1,z2+1.,n2mean-1.,beta2,2.0)
04452 + ecoul(z1+1.,n1mean-1.,beta1,z2-1.,n2mean+1.,beta2,2.0)
04453 - 2.0 * ecoul(z1,n1mean,beta1,z2,n2mean,beta2,2.0)
04454 - 2.0 * umass(z1, n1mean,beta1)
04455 - 2.0 * umass(z2, n2mean,beta2) ) * 0.5;
04456
04457
04458
04459
04460
04461
04462 za1width=std::sqrt( (0.5 * std::sqrt(1.0/a_levdens*(e_scission)) / cz) );
04463 za1width=std::sqrt( (max((za1width*za1width-(1.0/12.0)),0.1)) );
04464
04465
04466 a1r = z1 + n1mean;
04467 a1r = (float)((int)((a1r+0.5)));
04468 a2r = a - a1r;
04469
04470
04471
04472 n1ucd = n/a * a1r;
04473 n2ucd = n/a * a2r;
04474 z1ucd = z/a * a1r;
04475 z2ucd = z/a * a2r;
04476
04477 re1 = umass(z1ucd-1.,n1ucd+1.,beta1) + umass(z2ucd+1.,n2ucd-1.,beta2)
04478 + ecoul(z1ucd-1.,n1ucd+1.,beta1,z2ucd+1.,n2ucd-1.,beta2,d);
04479 re2 = umass(z1ucd,n1ucd,beta1) + umass(z2ucd,n2ucd,beta2)
04480 + ecoul(z1ucd,n1ucd,beta1,z2ucd,n2ucd,beta2,d);
04481 re3 = umass(z1ucd+1.,n1ucd-1.,beta1) + umass(z2ucd-1.,n2ucd+1.,beta2) +
04482 + ecoul(z1ucd+1.,n1ucd-1.,beta1,z2ucd-1.,n2ucd+1.,beta2,d);
04483
04484 eps2 = (re1-2.0*re2+re3) / 2.0;
04485 eps1 = (re3 - re1)/2.0;
04486 dn1_pol = - eps1 / (2.0 * eps2);
04487 z1 = z1ucd + dn1_pol;
04488 if (itest == 1) {
04489 G4cout << "before proton sample " << z1 << G4endl;
04490 }
04491
04492
04493 z1 = gausshaz(k, z1, za1width);
04494 if (itest == 1) {
04495 G4cout << "after proton sample " << z1 << G4endl;
04496 }
04497 even_odd(z1,r_e_o,i_help);
04498
04499 z1 = (float)(i_help);
04500 z2 = (float)((int)( (z - z1 + 0.5)) );
04501
04502 n1r = a1r - z1;
04503 n2r = n - n1r;
04504 }
04505 else {
04506
04507 cn = ( umass(z1, n1mean+1.,beta1) + umass(z1, n1mean-1., beta1)
04508 + umass(z2, n2mean+1.,beta2) + umass(z2, n2mean-1., beta2)
04509 + ecoul(z1,n1mean+1.,beta1,z2,n2mean-1.,beta2,2.0)
04510 + ecoul(z1,n1mean-1.,beta1,z2,n2mean+1.,beta2,2.0)
04511 - 2.0 * ecoul(z1,n1mean,beta1,z2,n2mean,beta2,2.0)
04512 - 2.0 * umass(z1, n1mean, 0.6)
04513 - 2.0 * umass(z2, n2mean, 0.6) ) * 0.5;
04514
04515
04516
04517
04518
04519
04520
04521 n1width=std::sqrt( (0.5 * std::sqrt(1.0/a_levdens*(e_scission)) / cn) );
04522 n1width=max(n1width, sigzmin);
04523
04524
04525
04526
04527 n1r = gausshaz(k, n1mean, n1width);
04528 n1r = (float)( (int)((n1r+0.5)) );
04529 n2r = n - n1r;
04530
04531 even_odd(z1,r_e_o,i_help);
04532
04533 z1 = (float)(i_help);
04534 z2 = z - z1;
04535
04536 a1r = z1 + n1r;
04537 a2r = z2 + n2r;
04538
04539 }
04540 }
04541
04542 if (itest == 1) {
04543 G4cout << "remid imode = " << imode << G4endl;
04544 G4cout << "n1width = " << n1width << G4endl;
04545 G4cout << "n1r = " << n1r << G4endl;
04546 G4cout << "a1r = " << a1r << G4endl;
04547 G4cout << "n2r = " << n2r << G4endl;
04548 G4cout << "a2r = " << a2r << G4endl;
04549 }
04550
04551
04552
04553 e1exc_sigma = 5.5;
04554 e2exc_sigma = 5.5;
04555
04556 neufcentquatrevingtsept:
04557
04558
04559
04560
04561 e1final = gausshaz(k, e1exc, e1exc_sigma);
04562 e2final = gausshaz(k, e2exc, e2exc_sigma);
04563 if ( (e1final < 0.0) || (e2final < 0.0) ) goto neufcentquatrevingtsept;
04564 if (itest == 1) {
04565 G4cout << "sampled exc 1 " << e1final << G4endl;
04566 G4cout << "sampled exc 2 " << e2final << G4endl;
04567 }
04568
04569
04570
04571
04572
04573
04574
04575
04576
04577
04578
04579
04580
04581
04582
04583
04584
04585
04586 n1 = n1r;
04587 n2 = n2r;
04588 a1 = n1 + z1;
04589 a2 = n2 + z2;
04590 e1 = e1final;
04591 e2 = e2final;
04592
04593
04594 tker = (z1 * z2 * 1.44) /
04595 ( r0 * std::pow(a1,0.33333) * (1.0 + 2.0/3.0 * beta1) +
04596 r0 * std::pow(a2,0.33333) * (1.0 + 2.0/3.0 * beta2) + 2.0 );
04597
04598 ekinr1 = tker * a2 / a;
04599
04600 ekinr2 = tker * a1 / a;
04601
04602 v1 = std::sqrt( (ekinr1/a1) ) * 1.3887;
04603 v2 = std::sqrt( (ekinr2/a2) ) * 1.3887;
04604
04605 if (itest == 1) {
04606 G4cout << "ekinr1 " << ekinr1 << G4endl;
04607 G4cout << "ekinr2 " << ekinr2 << G4endl;
04608 }
04609
04610 milledeux:
04611
04612
04613
04614
04615 if ( (icz == -1) || (a1 < 0.0) || (a2 < 0.0) ) {
04616
04617
04618
04619
04620
04621 if (itest == 1) {
04622 G4cout << "milledeux: liquid-drop option " << G4endl;
04623 }
04624
04625 n = a-z;
04626
04627 zsymm = z / 2.0;
04628 nsymm = n / 2.0;
04629 asymm = nsymm + zsymm;
04630
04631 a_levdens = a / xlevdens;
04632
04633 masscurv = 2.0;
04634 cz_symm = 8.0 / std::pow(z,2) * masscurv;
04635
04636 wzsymm = std::sqrt( (0.5 * std::sqrt(1.0/a_levdens*e) / cz_symm) ) ;
04637
04638 if (itest == 1) {
04639 G4cout << " symmetric high energy fission " << G4endl;
04640 G4cout << "wzsymm " << wzsymm << G4endl;
04641 }
04642
04643 z1mean = zsymm;
04644 z1width = wzsymm;
04645
04646
04647 z1 = 1.0;
04648 z2 = 1.0;
04649 while ( (z1 < 5.0) || (z2 < 5.0) ) {
04650
04651
04652 z1 = gausshaz(kkk, z1mean, z1width);
04653 z2 = z - z1;
04654 }
04655
04656 if (itest == 1) {
04657 G4cout << " z1 " << z1 << G4endl;
04658 G4cout << " z2 " << z2 << G4endl;
04659 }
04660 if (itest == 1) {
04661 G4cout << " zsymm " << zsymm << G4endl;
04662 G4cout << " nsymm " << nsymm << G4endl;
04663 G4cout << " asymm " << asymm << G4endl;
04664 }
04665
04666
04667
04668
04669
04670
04671
04672
04673 n1ucd = z1 * n/z;
04674 n2ucd = z2 * n/z;
04675 re1 = umass(z1,n1ucd,0.6) + umass(z2,n2ucd,0.6) +
04676 ecoul(z1,n1ucd,0.6,z2,n2ucd,0.6,2.0);
04677 re2 = umass(z1,n1ucd+1.,0.6) + umass(z2,n2ucd-1.,0.6) +
04678 ecoul(z1,n1ucd+1.,0.6,z2,n2ucd-1.,0.6,2.0);
04679 re3 = umass(z1,n1ucd+2.,0.6) + umass(z2,n2ucd-2.,0.6) +
04680 ecoul(z1,n1ucd+2.,0.6,z2,n2ucd-2.,0.6,2.0);
04681 reps2 = (re1-2.0*re2+re3)/2.0;
04682 reps1 = re2 - re1 -reps2;
04683 rn1_pol = -reps1/(2.0*reps2);
04684 n1mean = n1ucd + rn1_pol;
04685 n2mean = n - n1mean;
04686
04687 if (itest == 1) {
04688 G4cout << " n1mean " << n1mean << G4endl;
04689 G4cout << " n2mean " << n2mean << G4endl;
04690 }
04691
04692 cn = (umass(z1,n1mean+1.,0.0) + umass(z1,n1mean-1.,0.0) +
04693 + umass(z2,n2mean+1.,0.0) + umass(z2,n2mean-1.,0.0)
04694 - 2.0 * umass(z1,n1mean,0.0) +
04695 - 2.0 * umass(z2,n2mean,0.0) ) * 0.5;
04696
04697
04698 n1width = std::sqrt( (0.5 * std::sqrt(1.0/a_levdens*e) / cn) );
04699
04700 if (itest == 1) {
04701 G4cout << " cn " << cn << G4endl;
04702 G4cout << " n1width " << n1width << G4endl;
04703 }
04704
04705
04706
04707
04708 n1r = (float)( (int)(gausshaz(k, n1mean,n1width)) );
04709 n2r = n - n1r;
04710
04711 a1 = z1 + n1r;
04712 a2 = z2 + n2r;
04713
04714 e1 = e*a1/(a1+a2);
04715 e2 = e - e*a1/(a1+a2);
04716 if (itest == 1) {
04717 G4cout << " n1r " << n1r << G4endl;
04718 G4cout << " n2r " << n2r << G4endl;
04719 }
04720
04721 }
04722
04723 if (itest == 1) {
04724 G4cout << " a1 " << a1 << G4endl;
04725 G4cout << " z1 " << z1 << G4endl;
04726 G4cout << " a2 " << a2 << G4endl;
04727 G4cout << " z2 " << z2 << G4endl;
04728 G4cout << " e1 " << e1 << G4endl;
04729 G4cout << " e2 " << e << G4endl;
04730 }
04731
04732
04733 tker = (z1 * z2 * 1.44) /
04734 ( r0 * std::pow(a1,0.33333) * (1.0 + 2.0/3.0 * beta1) +
04735 r0 * std::pow(a2,0.33333) * (1.0 + 2.0/3.0 * beta2) + 2.0 );
04736
04737 ekin1 = tker * a2 / a;
04738
04739 ekin2 = tker * a1 / a;
04740
04741 v1 = std::sqrt( (ekin1/a1) ) * 1.3887;
04742 v2 = std::sqrt( (ekin2/a2) ) * 1.3887;
04743
04744 if (itest == 1) {
04745 G4cout << " kinetic energies " << G4endl;
04746 G4cout << " ekin1 " << ekin1 << G4endl;
04747 G4cout << " ekin2 " << ekin2 << G4endl;
04748 }
04749 }
04750
04751
04752 void G4Abla::translab(G4double gamrem, G4double etrem, G4double csrem[4], G4int nopart, G4int ndec)
04753 {
04754
04755
04756
04757
04758
04759
04760
04761
04762
04763
04764
04765
04766
04767
04768
04769
04770
04771
04772
04773
04774
04775
04776
04777
04778
04779
04780
04781 G4double sitet = std::sqrt(std::pow(csrem[1],2)+std::pow(csrem[2],2));
04782 G4double cstet = 0.0, siphi = 0.0, csphi = 0.0;
04783 G4double R[4][4];
04784 for(G4int init_i = 0; init_i < 4; init_i++) {
04785 for(G4int init_j = 0; init_j < 4; init_j++) {
04786 R[init_i][init_j] = 0.0;
04787 }
04788 }
04789
04790 if(sitet > 1.0e-6) {
04791 cstet = csrem[3];
04792 siphi = csrem[2]/sitet;
04793 csphi = csrem[1]/sitet;
04794
04795 R[1][1] = cstet*csphi;
04796 R[1][2] = -siphi;
04797 R[1][3] = sitet*csphi;
04798 R[2][1] = cstet*siphi;
04799 R[2][2] = csphi;
04800 R[2][3] = sitet*siphi;
04801 R[3][1] = -sitet;
04802 R[3][2] = 0.0;
04803 R[3][3] = cstet;
04804 }
04805 else {
04806 R[1][1] = 1.0;
04807 R[1][2] = 0.0;
04808 R[1][3] = 0.0;
04809 R[2][1] = 0.0;
04810 R[2][2] = 1.0;
04811 R[2][3] = 0.0;
04812 R[3][1] = 0.0;
04813 R[3][2] = 0.0;
04814 R[3][3] = 1.0;
04815 }
04816
04817 G4int intp = 0;
04818 G4double el = 0.0;
04819 G4double masse = 0.0;
04820 G4double er = 0.0;
04821 G4double plabi[4];
04822 G4double ptrav2 = 0.0;
04823 G4double plabf[4];
04824 G4double bidon = 0.0;
04825 for(G4int init_i = 0; init_i < 4; init_i++) {
04826 plabi[init_i] = 0.0;
04827 plabf[init_i] = 0.0;
04828 }
04829
04830 for(G4int i = ndec; i <= volant->iv; i++) {
04831 intp = i + nopart;
04832 varntp->ntrack = varntp->ntrack + 1;
04833 if(nint(volant->acv[i]) == 0 && nint(volant->zpcv[i]) == 0) {
04834 if(verboseLevel > 2) {
04835 G4cout <<"Error: Particles with A = 0 Z = 0 detected! " << G4endl;
04836 }
04837 continue;
04838 }
04839 if(varntp->ntrack >= VARNTPSIZE) {
04840 if(verboseLevel > 2) {
04841 G4cout <<"Error! Output data structure not big enough!" << G4endl;
04842 }
04843 }
04844 varntp->avv[intp] = nint(volant->acv[i]);
04845 varntp->zvv[intp] = nint(volant->zpcv[i]);
04846 varntp->itypcasc[intp] = 0;
04847
04848 if (varntp->avv[intp] == -1) {
04849 masse = 138.00;
04850
04851
04852 }
04853 else {
04854 mglms(double(volant->acv[i]),double(volant->zpcv[i]),0, &el);
04855
04856 masse = volant->zpcv[i]*938.27 + (volant->acv[i] - volant->zpcv[i])*939.56 + el;
04857 }
04858
04859 er = std::sqrt(std::pow(volant->pcv[i],2) + std::pow(masse,2));
04860
04861 plabi[1] = volant->pcv[i]*(volant->xcv[i]);
04862 plabi[2] = volant->pcv[i]*(volant->ycv[i]);
04863 plabi[3] = er*etrem + gamrem*(volant->pcv[i])*(volant->zcv[i]);
04864
04865 ptrav2 = std::pow(plabi[1],2) + std::pow(plabi[2],2) + std::pow(plabi[3],2);
04866
04867 varntp->plab[intp] = std::sqrt(ptrav2);
04868 varntp->enerj[intp] = std::sqrt(ptrav2 + std::pow(masse,2)) - masse;
04869
04870
04871 for(G4int j = 1; j <= 3; j++) {
04872 plabf[j] = 0.0;
04873 for(G4int k = 1; k <= 3; k++) {
04874 plabf[j] = plabf[j] + R[k][j]*plabi[k];
04875 }
04876 }
04877
04878 volant->pcv[i] = varntp->plab[intp];
04879 ptrav2 = std::sqrt(std::pow(plabf[1],2) + std::pow(plabf[2],2) + std::pow(plabf[3],2));
04880 if(ptrav2 >= 1.0e-6) {
04881 volant->xcv[i] = plabf[1]/ptrav2;
04882 volant->ycv[i] = plabf[2]/ptrav2;
04883 volant->zcv[i] = plabf[3]/ptrav2;
04884 }
04885 else {
04886 volant->xcv[i] = 1.0;
04887 volant->ycv[i] = 0.0;
04888 volant->zcv[i] = 0.0;
04889 }
04890
04891 if(varntp->plab[intp] >= 1.0e-6) {
04892 bidon = plabf[3]/(varntp->plab[intp]);
04893
04894 if(bidon > 1.0) {
04895 bidon = 1.0;
04896 }
04897 if(bidon < -1.0) {
04898 bidon = -1.0;
04899 }
04900 varntp->tetlab[intp] = std::acos(bidon);
04901 sitet = std::sin(varntp->tetlab[intp]);
04902 varntp->philab[intp] = std::atan2(plabf[2],plabf[1]);
04903 varntp->tetlab[intp] = varntp->tetlab[intp]*57.2957795;
04904 varntp->philab[intp] = varntp->philab[intp]*57.2957795;
04905 }
04906 else {
04907 varntp->tetlab[intp] = 90.0;
04908 varntp->philab[intp] = 0.0;
04909 }
04910 }
04911 }
04912
04913
04914
04915
04916 void G4Abla::translabpf(G4double masse1, G4double t1, G4double p1, G4double ctet1,
04917 G4double phi1, G4double gamrem, G4double etrem, G4double R[][4],
04918 G4double *plab1, G4double *gam1, G4double *eta1, G4double csdir[])
04919 {
04920
04921
04922
04923
04924
04925
04926
04927
04928
04929
04930
04931
04932
04933 G4double er = t1 + masse1;
04934
04935 G4double sitet = std::sqrt(1.0 - std::pow(ctet1,2));
04936
04937 G4double plabi[4];
04938 G4double plabf[4];
04939 for(G4int init_i = 0; init_i < 4; init_i++) {
04940 plabi[init_i] = 0.0;
04941 plabf[init_i] = 0.0;
04942 }
04943
04944
04945 plabi[1] = p1*sitet*std::cos(phi1);
04946 plabi[2] = p1*sitet*std::sin(phi1);
04947 plabi[3] = er*etrem + gamrem*p1*ctet1;
04948
04949
04950 for(G4int j = 1; j <= 3; j++) {
04951 plabf[j] = 0.0;
04952 for(G4int k = 1; k <= 3; k++) {
04953
04954 plabf[j] = plabf[j] + R[k][j]*plabi[k];
04955 }
04956 }
04957
04958
04959 (*plab1) = std::pow(plabf[1],2) + std::pow(plabf[2],2) + std::pow(plabf[3],2);
04960 (*gam1) = std::sqrt(std::pow(masse1,2) + (*plab1))/masse1;
04961 (*plab1) = std::sqrt((*plab1));
04962 (*eta1) = (*plab1)/masse1;
04963
04964 if((*plab1) <= 1.0e-6) {
04965 csdir[1] = 0.0;
04966 csdir[2] = 0.0;
04967 csdir[3] = 1.0;
04968 }
04969 else {
04970 for(G4int i = 1; i <= 3; i++) {
04971 csdir[i] = plabf[i]/(*plab1);
04972 }
04973 }
04974 }
04975
04976
04977 void G4Abla::lorab(G4double gam, G4double eta, G4double ein, G4double pin[],
04978 G4double *eout, G4double pout[])
04979 {
04980
04981
04982
04983
04984
04985
04986 pout[1] = pin[1];
04987 pout[2] = pin[2];
04988 (*eout) = gam*ein + eta*pin[3];
04989 pout[3] = eta*ein + gam*pin[3];
04990 }
04991
04992
04993 void G4Abla::rotab(G4double R[4][4], G4double pin[4], G4double pout[4])
04994 {
04995
04996
04997
04998
04999 for(G4int i = 1; i <= 3; i++) {
05000 pout[i] = 0.0;
05001 for(G4int j = 1; j <= 3; j++) {
05002
05003 pout[i] = pout[i] + R[j][i]*pin[j];
05004 }
05005 }
05006 }
05007
05008
05009
05010
05011
05012 void G4Abla::standardRandom(G4double *rndm, G4long *seed)
05013 {
05014 (*seed) = (*seed);
05015
05016 (*rndm) = G4UniformRand();
05017 }
05018
05019 G4double G4Abla::haz(G4int k)
05020 {
05021 const G4int pSize = 110;
05022 static G4double p[pSize];
05023 static G4long ix = 0, i = 0;
05024 static G4double x = 0.0, y = 0.0, a = 0.0, haz = 0.0;
05025
05026
05027
05028
05029
05030 if(ix == 0) {
05031 ix = hazard->ial;
05032 }
05033
05034 if (k <= -1) {
05035 if(k == -1) {
05036 ix = 0;
05037 }
05038 else {
05039 x = 0.0;
05040 y = secnds(int(x));
05041 ix = int(y * 100 + 43543000);
05042 if(mod(ix,2) == 0) {
05043 ix = ix + 1;
05044 }
05045 }
05046
05047
05048
05049
05050
05051
05052 standardRandom(&x, &ix);
05053 for(G4int i = 0; i < pSize; i++) {
05054 standardRandom(&(p[i]), &ix);
05055 }
05056 standardRandom(&a, &ix);
05057 k = 0;
05058 }
05059
05060 i = nint(100*a)+1;
05061 haz = p[i];
05062 standardRandom(&a, &ix);
05063 p[i] = a;
05064
05065 hazard->ial = ix;
05066
05067 return haz;
05068 }
05069
05070
05071 G4double G4Abla::gausshaz(int k, double xmoy, double sig)
05072 {
05073
05074
05075
05076 static G4int iset = 0;
05077 static G4double v1,v2,r,fac,gset,gausshaz;
05078
05079 if(iset == 0) {
05080 do {
05081 v1 = 2.0*haz(k) - 1.0;
05082 v2 = 2.0*haz(k) - 1.0;
05083 r = std::pow(v1,2) + std::pow(v2,2);
05084 } while(r >= 1);
05085
05086 fac = std::sqrt(-2.*std::log(r)/r);
05087
05088 gset = v1*fac;
05089 gausshaz = v2*fac*sig+xmoy;
05090 iset = 1;
05091 }
05092 else {
05093 gausshaz=gset*sig+xmoy;
05094 iset=0;
05095 }
05096
05097 return gausshaz;
05098 }
05099
05100
05101
05102
05103 G4double G4Abla::min(G4double a, G4double b)
05104 {
05105 if(a < b) {
05106 return a;
05107 }
05108 else {
05109 return b;
05110 }
05111 }
05112
05113 G4int G4Abla::min(G4int a, G4int b)
05114 {
05115 if(a < b) {
05116 return a;
05117 }
05118 else {
05119 return b;
05120 }
05121 }
05122
05123 G4double G4Abla::max(G4double a, G4double b)
05124 {
05125 if(a > b) {
05126 return a;
05127 }
05128 else {
05129 return b;
05130 }
05131 }
05132
05133 G4int G4Abla::max(G4int a, G4int b)
05134 {
05135 if(a > b) {
05136 return a;
05137 }
05138 else {
05139 return b;
05140 }
05141 }
05142
05143 G4int G4Abla::nint(G4double number)
05144 {
05145 G4double intpart = 0.0;
05146 G4double fractpart = 0.0;
05147 fractpart = std::modf(number, &intpart);
05148 if(number == 0) {
05149 return 0;
05150 }
05151 if(number > 0) {
05152 if(fractpart < 0.5) {
05153 return int(std::floor(number));
05154 }
05155 else {
05156 return int(std::ceil(number));
05157 }
05158 }
05159 if(number < 0) {
05160 if(fractpart < -0.5) {
05161 return int(std::floor(number));
05162 }
05163 else {
05164 return int(std::ceil(number));
05165 }
05166 }
05167
05168 return int(std::floor(number));
05169 }
05170
05171 G4int G4Abla::secnds(G4int x)
05172 {
05173 time_t mytime;
05174 tm *mylocaltime;
05175
05176 time(&mytime);
05177 mylocaltime = localtime(&mytime);
05178
05179 if(x == 0) {
05180 return(mylocaltime->tm_hour*60*60 + mylocaltime->tm_min*60 + mylocaltime->tm_sec);
05181 }
05182 else {
05183 return(mytime - x);
05184 }
05185 }
05186
05187 G4int G4Abla::mod(G4int a, G4int b)
05188 {
05189 if(b != 0) {
05190 return (a - (a/b)*b);
05191 }
05192 else {
05193 return 0;
05194 }
05195 }
05196
05197 G4double G4Abla::dmod(G4double a, G4double b)
05198 {
05199 if(b != 0) {
05200 return (a - (a/b)*b);
05201 }
05202 else {
05203 return 0.0;
05204 }
05205 }
05206
05207 G4double G4Abla::dint(G4double a)
05208 {
05209 G4double value = 0.0;
05210
05211 if(a < 0.0) {
05212 value = double(std::ceil(a));
05213 }
05214 else {
05215 value = double(std::floor(a));
05216 }
05217
05218 return value;
05219 }
05220
05221 G4int G4Abla::idint(G4double a)
05222 {
05223 G4int value = 0;
05224
05225 if(a < 0) {
05226 value = int(std::ceil(a));
05227 }
05228 else {
05229 value = int(std::floor(a));
05230 }
05231
05232 return value;
05233 }
05234
05235 G4int G4Abla::idnint(G4double value)
05236 {
05237 G4double valueCeil = int(std::ceil(value));
05238 G4double valueFloor = int(std::floor(value));
05239
05240 if(std::fabs(value - valueCeil) < std::fabs(value - valueFloor)) {
05241 return int(valueCeil);
05242 }
05243 else {
05244 return int(valueFloor);
05245 }
05246 }
05247
05248 G4double G4Abla::dmin1(G4double a, G4double b, G4double c)
05249 {
05250 if(a < b && a < c) {
05251 return a;
05252 }
05253 if(b < a && b < c) {
05254 return b;
05255 }
05256 if(c < a && c < b) {
05257 return c;
05258 }
05259 return a;
05260 }
05261
05262 G4double G4Abla::utilabs(G4double a)
05263 {
05264 if(a > 0) {
05265 return a;
05266 }
05267 if(a < 0) {
05268 return (-1*a);
05269 }
05270 if(a == 0) {
05271 return a;
05272 }
05273
05274 return a;
05275 }