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 #include <string.h>
00031 #include <cmath>
00032 #include <stdio.h>
00033
00034 #include "G4AntiNeutronAnnihilationAtRest.hh"
00035 #include "G4SystemOfUnits.hh"
00036 #include "G4DynamicParticle.hh"
00037 #include "G4ParticleTypes.hh"
00038 #include "G4HadronicProcessStore.hh"
00039 #include "G4HadronicDeprecate.hh"
00040 #include "Randomize.hh"
00041
00042 #define MAX_SECONDARIES 100
00043
00044
00045
00046 G4AntiNeutronAnnihilationAtRest::G4AntiNeutronAnnihilationAtRest(const G4String& processName,
00047 G4ProcessType aType) :
00048 G4VRestProcess (processName, aType),
00049 massPionMinus(G4PionMinus::PionMinus()->GetPDGMass()/GeV),
00050 massPionZero(G4PionZero::PionZero()->GetPDGMass()/GeV),
00051 massPionPlus(G4PionPlus::PionPlus()->GetPDGMass()/GeV),
00052 massGamma(G4Gamma::Gamma()->GetPDGMass()/GeV),
00053 massAntiNeutron(G4AntiNeutron::AntiNeutron()->GetPDGMass()/GeV),
00054 massNeutron(G4Neutron::Neutron()->GetPDGMass()/GeV),
00055 pdefGamma(G4Gamma::Gamma()),
00056 pdefPionPlus(G4PionPlus::PionPlus()),
00057 pdefPionZero(G4PionZero::PionZero()),
00058 pdefPionMinus(G4PionMinus::PionMinus()),
00059 pdefProton(G4Proton::Proton()),
00060 pdefNeutron(G4Neutron::Neutron()),
00061 pdefAntiNeutron(G4AntiNeutron::AntiNeutron()),
00062 pdefDeuteron(G4Deuteron::Deuteron()),
00063 pdefTriton(G4Triton::Triton()),
00064 pdefAlpha(G4Alpha::Alpha())
00065 {
00066 G4HadronicDeprecate("G4AntiNeutronAnnihilationAtRest");
00067 if (verboseLevel>0) {
00068 G4cout << GetProcessName() << " is created "<< G4endl;
00069 }
00070 SetProcessSubType(fHadronAtRest);
00071 pv = new G4GHEKinematicsVector [MAX_SECONDARIES+1];
00072 eve = new G4GHEKinematicsVector [MAX_SECONDARIES];
00073 gkin = new G4GHEKinematicsVector [MAX_SECONDARIES];
00074
00075 G4HadronicProcessStore::Instance()->RegisterExtraProcess(this);
00076 }
00077
00078
00079
00080 G4AntiNeutronAnnihilationAtRest::~G4AntiNeutronAnnihilationAtRest()
00081 {
00082 G4HadronicProcessStore::Instance()->DeRegisterExtraProcess(this);
00083 delete [] pv;
00084 delete [] eve;
00085 delete [] gkin;
00086 }
00087
00088 void G4AntiNeutronAnnihilationAtRest::PreparePhysicsTable(const G4ParticleDefinition& p)
00089 {
00090 G4HadronicProcessStore::Instance()->RegisterParticleForExtraProcess(this, &p);
00091 }
00092
00093 void G4AntiNeutronAnnihilationAtRest::BuildPhysicsTable(const G4ParticleDefinition& p)
00094 {
00095 G4HadronicProcessStore::Instance()->PrintInfo(&p);
00096 }
00097
00098
00099
00100 G4bool G4AntiNeutronAnnihilationAtRest::IsApplicable(
00101 const G4ParticleDefinition& particle
00102 )
00103 {
00104 return ( &particle == pdefAntiNeutron );
00105
00106 }
00107
00108
00109 G4int G4AntiNeutronAnnihilationAtRest::GetNumberOfSecondaries()
00110 {
00111 return ( ngkine );
00112
00113 }
00114
00115
00116 G4GHEKinematicsVector* G4AntiNeutronAnnihilationAtRest::GetSecondaryKinematics()
00117 {
00118 return ( &gkin[0] );
00119
00120 }
00121
00122 G4double G4AntiNeutronAnnihilationAtRest::AtRestGetPhysicalInteractionLength(
00123 const G4Track& track,
00124 G4ForceCondition* condition
00125 )
00126 {
00127
00128 ResetNumberOfInteractionLengthLeft();
00129
00130
00131 *condition = NotForced;
00132
00133
00134 currentInteractionLength = GetMeanLifeTime(track, condition);
00135
00136 if ((currentInteractionLength <0.0) || (verboseLevel>2)){
00137 G4cout << "G4AntiNeutronAnnihilationAtRestProcess::AtRestGetPhysicalInteractionLength ";
00138 G4cout << "[ " << GetProcessName() << "]" <<G4endl;
00139 track.GetDynamicParticle()->DumpInfo();
00140 G4cout << " in Material " << track.GetMaterial()->GetName() <<G4endl;
00141 G4cout << "MeanLifeTime = " << currentInteractionLength/ns << "[ns]" <<G4endl;
00142 }
00143
00144 return theNumberOfInteractionLengthLeft * currentInteractionLength;
00145
00146 }
00147
00148 G4VParticleChange* G4AntiNeutronAnnihilationAtRest::AtRestDoIt(
00149 const G4Track& track,
00150 const G4Step&
00151 )
00152
00153
00154
00155
00156
00157 {
00158
00159
00160
00161
00162
00163 aParticleChange.Initialize(track);
00164
00165
00166
00167 globalTime = track.GetGlobalTime()/s;
00168 G4Material * aMaterial = track.GetMaterial();
00169 const G4int numberOfElements = aMaterial->GetNumberOfElements();
00170 const G4ElementVector* theElementVector = aMaterial->GetElementVector();
00171
00172 const G4double* theAtomicNumberDensity = aMaterial->GetAtomicNumDensityVector();
00173 G4double normalization = 0;
00174 for ( G4int i1=0; i1 < numberOfElements; i1++ )
00175 {
00176 normalization += theAtomicNumberDensity[i1] ;
00177
00178 }
00179 G4double runningSum= 0.;
00180 G4double random = G4UniformRand()*normalization;
00181 for ( G4int i2=0; i2 < numberOfElements; i2++ )
00182 {
00183 runningSum += theAtomicNumberDensity[i2];
00184
00185 if (random<=runningSum)
00186 {
00187 targetCharge = G4double( ((*theElementVector)[i2])->GetZ());
00188 targetAtomicMass = (*theElementVector)[i2]->GetN();
00189 }
00190 }
00191 if (random>runningSum)
00192 {
00193 targetCharge = G4double( ((*theElementVector)[numberOfElements-1])->GetZ());
00194 targetAtomicMass = (*theElementVector)[numberOfElements-1]->GetN();
00195 }
00196
00197 if (verboseLevel>1) {
00198 G4cout << "G4AntiNeutronAnnihilationAtRest::AtRestDoIt is invoked " <<G4endl;
00199 }
00200
00201 G4ParticleMomentum momentum;
00202 G4float localtime;
00203
00204 G4ThreeVector position = track.GetPosition();
00205
00206 GenerateSecondaries();
00207
00208 aParticleChange.SetNumberOfSecondaries( ngkine );
00209
00210 for ( G4int isec = 0; isec < ngkine; isec++ ) {
00211 G4DynamicParticle* aNewParticle = new G4DynamicParticle;
00212 aNewParticle->SetDefinition( gkin[isec].GetParticleDef() );
00213 aNewParticle->SetMomentum( gkin[isec].GetMomentum() * GeV );
00214
00215 localtime = globalTime + gkin[isec].GetTOF();
00216
00217 G4Track* aNewTrack = new G4Track( aNewParticle, localtime*s, position );
00218 aNewTrack->SetTouchableHandle(track.GetTouchableHandle());
00219 aParticleChange.AddSecondary( aNewTrack );
00220
00221 }
00222
00223 aParticleChange.ProposeLocalEnergyDeposit( 0.0*GeV );
00224
00225 aParticleChange.ProposeTrackStatus(fStopAndKill);
00226
00227
00228
00229 ResetNumberOfInteractionLengthLeft();
00230
00231 return &aParticleChange;
00232
00233 }
00234
00235
00236 void G4AntiNeutronAnnihilationAtRest::GenerateSecondaries()
00237 {
00238 static G4int index;
00239 static G4int l;
00240 static G4int nopt;
00241 static G4int i;
00242
00243
00244 for (i = 1; i <= MAX_SECONDARIES; ++i) {
00245 pv[i].SetZero();
00246 }
00247
00248
00249 ngkine = 0;
00250 ntot = 0;
00251 result.SetZero();
00252 result.SetMass( massAntiNeutron );
00253 result.SetKineticEnergyAndUpdate( 0. );
00254 result.SetTOF( 0. );
00255 result.SetParticleDef( pdefAntiNeutron );
00256
00257
00258
00259 AntiNeutronAnnihilation(&nopt);
00260
00261
00262 if (ntot != 0 || result.GetParticleDef() != pdefAntiNeutron) {
00263
00264
00265
00266
00267
00268
00269
00270 gkin[0] = result;
00271 gkin[0].SetTOF( result.GetTOF() * 5e-11 );
00272 ngkine = 1;
00273
00274
00275
00276
00277
00278 for (l = 1; l <= ntot; ++l) {
00279 index = l - 1;
00280
00281
00282
00283 if (ngkine < MAX_SECONDARIES) {
00284 gkin[ngkine] = eve[index];
00285 gkin[ngkine].SetTOF( eve[index].GetTOF() * 5e-11 );
00286 ++ngkine;
00287 }
00288 }
00289 }
00290 else {
00291
00292
00293 ngkine = 0;
00294 ntot = 0;
00295 globalTime += result.GetTOF() * G4float(5e-11);
00296 }
00297
00298
00299 ngkine = G4int(std::min(ngkine,G4int(MAX_SECONDARIES)));
00300
00301 }
00302
00303
00304 void G4AntiNeutronAnnihilationAtRest::Poisso(G4float xav, G4int *iran)
00305 {
00306 static G4int i;
00307 static G4float r, p1, p2, p3;
00308 static G4int fivex;
00309 static G4float rr, ran, rrr, ran1;
00310
00311
00312
00313
00314
00315
00316 if (xav > G4float(9.9)) {
00317
00318 Normal(&ran1);
00319 ran1 = xav + ran1 * std::sqrt(xav);
00320 *iran = G4int(ran1);
00321 if (*iran < 0) {
00322 *iran = 0;
00323 }
00324 }
00325 else {
00326 fivex = G4int(xav * G4float(5.));
00327 *iran = 0;
00328 if (fivex > 0) {
00329 r = std::exp(-G4double(xav));
00330 ran1 = G4UniformRand();
00331 if (ran1 > r) {
00332 rr = r;
00333 for (i = 1; i <= fivex; ++i) {
00334 ++(*iran);
00335 if (i <= 5) {
00336 rrr = std::pow(xav, G4float(i)) / NFac(i);
00337 }
00338
00339 if (i > 5) {
00340 rrr = std::exp(i * std::log(xav) -
00341 (i + G4float(.5)) * std::log(i * G4float(1.)) +
00342 i - G4float(.9189385));
00343 }
00344 rr += r * rrr;
00345 if (ran1 <= rr) {
00346 break;
00347 }
00348 }
00349 }
00350 }
00351 else {
00352
00353 p1 = xav * std::exp(-G4double(xav));
00354 p2 = xav * p1 / G4float(2.);
00355 p3 = xav * p2 / G4float(3.);
00356 ran = G4UniformRand();
00357 if (ran >= p3) {
00358 if (ran >= p2) {
00359 if (ran >= p1) {
00360 *iran = 0;
00361 }
00362 else {
00363 *iran = 1;
00364 }
00365 }
00366 else {
00367 *iran = 2;
00368 }
00369 }
00370 else {
00371 *iran = 3;
00372 }
00373 }
00374 }
00375
00376 }
00377
00378
00379 G4int G4AntiNeutronAnnihilationAtRest::NFac(G4int n)
00380 {
00381 G4int ret_val;
00382
00383 static G4int i, j;
00384
00385
00386
00387
00388 ret_val = 1;
00389 j = n;
00390 if (j > 1) {
00391 if (j > 10) {
00392 j = 10;
00393 }
00394 for (i = 2; i <= j; ++i) {
00395 ret_val *= i;
00396 }
00397 }
00398 return ret_val;
00399
00400 }
00401
00402
00403 void G4AntiNeutronAnnihilationAtRest::Normal(G4float *ran)
00404 {
00405 static G4int i;
00406
00407
00408
00409
00410 *ran = G4float(-6.);
00411 for (i = 1; i <= 12; ++i) {
00412 *ran += G4UniformRand();
00413 }
00414
00415 }
00416
00417
00418 void G4AntiNeutronAnnihilationAtRest::AntiNeutronAnnihilation(G4int *nopt)
00419 {
00420 static G4float brr[3] = { G4float(.125),G4float(.25),G4float(.5) };
00421
00422 G4float r__1;
00423
00424 static G4int i, ii, kk;
00425 static G4int nt;
00426 static G4float cfa, eka;
00427 static G4int ika, nbl;
00428 static G4float ran, pcm;
00429 static G4int isw;
00430 static G4float tex;
00431 static G4ParticleDefinition* ipa1;
00432 static G4float ran1, ran2, ekin, tkin;
00433 static G4float targ;
00434 static G4ParticleDefinition* inve;
00435 static G4float ekin1, ekin2, black;
00436 static G4float pnrat, rmnve1, rmnve2;
00437 static G4float ek, en;
00438
00439
00440
00441
00442
00443
00444
00445
00446
00447
00448
00449 pv[1].SetZero();
00450 pv[1].SetMass( massAntiNeutron );
00451 pv[1].SetKineticEnergyAndUpdate( 0. );
00452 pv[1].SetTOF( result.GetTOF() );
00453 pv[1].SetParticleDef( result.GetParticleDef() );
00454 isw = 1;
00455 ran = G4UniformRand();
00456 if (ran > brr[0]) {
00457 isw = 2;
00458 }
00459 if (ran > brr[1]) {
00460 isw = 3;
00461 }
00462 if (ran > brr[2]) {
00463 isw = 4;
00464 }
00465 *nopt = isw;
00466
00467
00468
00469 rmnve1 = massPionPlus;
00470 rmnve2 = massPionMinus;
00471 if (isw == 2) {
00472 rmnve1 = massPionZero;
00473 }
00474 if (isw == 2) {
00475 rmnve2 = massPionZero;
00476 }
00477 if (isw == 3) {
00478 rmnve2 = massPionZero;
00479 }
00480 if (isw == 4) {
00481 rmnve1 = massGamma;
00482 }
00483 if (isw == 4) {
00484 rmnve2 = massGamma;
00485 }
00486 ek = massNeutron + massAntiNeutron - rmnve1 - rmnve2;
00487 tkin = ExNu(ek);
00488 ek -= tkin;
00489 if (ek < G4float(1e-4)) {
00490 ek = G4float(1e-4);
00491 }
00492 ek /= G4float(2.);
00493 en = ek + (rmnve1 + rmnve2) / G4float(2.);
00494 r__1 = en * en - rmnve1 * rmnve2;
00495 pcm = r__1 > 0 ? std::sqrt(r__1) : 0;
00496 pv[2].SetZero();
00497 pv[2].SetMass( rmnve1 );
00498 pv[3].SetZero();
00499 pv[3].SetMass( rmnve2 );
00500 if (isw > 3) {
00501 pv[2].SetMass( 0. );
00502 pv[3].SetMass( 0. );
00503 }
00504 pv[2].SetEnergyAndUpdate( std::sqrt(pv[2].GetMass()*pv[2].GetMass()+pcm*pcm) );
00505 pv[2].SetTOF( result.GetTOF() );
00506 pv[3].SetEnergy( std::sqrt(pv[3].GetMass()*pv[3].GetMass()+pcm*pcm) );
00507 pv[3].SetMomentumAndUpdate( -pv[2].GetMomentum().x(), -pv[2].GetMomentum().y(), -pv[2].GetMomentum().z() );
00508 pv[3].SetTOF( result.GetTOF() );
00509 switch ((int)isw) {
00510 case 1:
00511 pv[2].SetParticleDef( pdefPionPlus );
00512 pv[3].SetParticleDef( pdefPionMinus );
00513 break;
00514 case 2:
00515 pv[2].SetParticleDef( pdefPionZero );
00516 pv[3].SetParticleDef( pdefPionZero );
00517 break;
00518 case 3:
00519 pv[2].SetParticleDef( pdefPionPlus );
00520 pv[3].SetParticleDef( pdefPionZero );
00521 break;
00522 case 4:
00523 pv[2].SetParticleDef( pdefGamma );
00524 pv[3].SetParticleDef( pdefGamma );
00525 break;
00526 default:
00527 break;
00528 }
00529 nt = 3;
00530 if (targetAtomicMass >= G4float(1.5)) {
00531 cfa = (targetAtomicMass - G4float(1.)) / G4float(120.) *
00532 G4float(.025) * std::exp(-G4double(targetAtomicMass - G4float(1.)) /
00533 G4float(120.));
00534 targ = G4float(1.);
00535 tex = evapEnergy1;
00536 if (tex >= G4float(.001)) {
00537 black = (targ * G4float(1.25) +
00538 G4float(1.5)) * evapEnergy1 / (evapEnergy1 + evapEnergy3);
00539 Poisso(black, &nbl);
00540 if (G4float(G4int(targ) + nbl) > targetAtomicMass) {
00541 nbl = G4int(targetAtomicMass - targ);
00542 }
00543 if (nt + nbl > (MAX_SECONDARIES - 2)) {
00544 nbl = (MAX_SECONDARIES - 2) - nt;
00545 }
00546 if (nbl > 0) {
00547 ekin = tex / nbl;
00548 ekin2 = G4float(0.);
00549 for (i = 1; i <= nbl; ++i) {
00550 if (nt == (MAX_SECONDARIES - 2)) {
00551 continue;
00552 }
00553 if (ekin2 > tex) {
00554 break;
00555 }
00556 ran1 = G4UniformRand();
00557 Normal(&ran2);
00558 ekin1 = -G4double(ekin) * std::log(ran1) -
00559 cfa * (ran2 * G4float(.5) + G4float(1.));
00560 if (ekin1 < G4float(0.)) {
00561 ekin1 = std::log(ran1) * G4float(-.01);
00562 }
00563 ekin1 *= G4float(1.);
00564 ekin2 += ekin1;
00565 if (ekin2 > tex) {
00566 ekin1 = tex - (ekin2 - ekin1);
00567 }
00568 if (ekin1 < G4float(0.)) {
00569 ekin1 = G4float(.001);
00570 }
00571 ipa1 = pdefNeutron;
00572 pnrat = G4float(1.) - targetCharge / targetAtomicMass;
00573 if (G4UniformRand() > pnrat) {
00574 ipa1 = pdefProton;
00575 }
00576 ++nt;
00577 pv[nt].SetZero();
00578 pv[nt].SetMass( ipa1->GetPDGMass()/GeV );
00579 pv[nt].SetKineticEnergyAndUpdate( ekin1 );
00580 pv[nt].SetTOF( result.GetTOF() );
00581 pv[nt].SetParticleDef( ipa1 );
00582 }
00583 if (targetAtomicMass >= G4float(230.) && ek <= G4float(2.)) {
00584 ii = nt + 1;
00585 kk = 0;
00586 eka = ek;
00587 if (eka > G4float(1.)) {
00588 eka *= eka;
00589 }
00590 if (eka < G4float(.1)) {
00591 eka = G4float(.1);
00592 }
00593 ika = G4int(G4float(3.6) / eka);
00594 for (i = 1; i <= nt; ++i) {
00595 --ii;
00596 if (pv[ii].GetParticleDef() != pdefProton) {
00597 continue;
00598 }
00599 ipa1 = pdefNeutron;
00600 pv[ii].SetMass( ipa1->GetPDGMass()/GeV );
00601 pv[ii].SetParticleDef( ipa1 );
00602 ++kk;
00603 if (kk > ika) {
00604 break;
00605 }
00606 }
00607 }
00608 }
00609 }
00610
00611
00612
00613 tex = evapEnergy3;
00614 if (tex >= G4float(.001)) {
00615 black = (targ * G4float(1.25) + G4float(1.5)) * evapEnergy3 /
00616 (evapEnergy1 + evapEnergy3);
00617 Poisso(black, &nbl);
00618 if (nt + nbl > (MAX_SECONDARIES - 2)) {
00619 nbl = (MAX_SECONDARIES - 2) - nt;
00620 }
00621 if (nbl > 0) {
00622 ekin = tex / nbl;
00623 ekin2 = G4float(0.);
00624 for (i = 1; i <= nbl; ++i) {
00625 if (nt == (MAX_SECONDARIES - 2)) {
00626 continue;
00627 }
00628 if (ekin2 > tex) {
00629 break;
00630 }
00631 ran1 = G4UniformRand();
00632 Normal(&ran2);
00633 ekin1 = -G4double(ekin) * std::log(ran1) -
00634 cfa * (ran2 * G4float(.5) + G4float(1.));
00635 if (ekin1 < G4float(0.)) {
00636 ekin1 = std::log(ran1) * G4float(-.01);
00637 }
00638 ekin1 *= G4float(1.);
00639 ekin2 += ekin1;
00640 if (ekin2 > tex) {
00641 ekin1 = tex - (ekin2 - ekin1);
00642 }
00643 if (ekin1 < G4float(0.)) {
00644 ekin1 = G4float(.001);
00645 }
00646 ran = G4UniformRand();
00647 inve = pdefDeuteron;
00648 if (ran > G4float(.6)) {
00649 inve = pdefTriton;
00650 }
00651 if (ran > G4float(.9)) {
00652 inve = pdefAlpha;
00653 }
00654 ++nt;
00655 pv[nt].SetZero();
00656 pv[nt].SetMass( inve->GetPDGMass()/GeV );
00657 pv[nt].SetKineticEnergyAndUpdate( ekin1 );
00658 pv[nt].SetTOF( result.GetTOF() );
00659 pv[nt].SetParticleDef( inve );
00660 }
00661 }
00662 }
00663 }
00664 result = pv[2];
00665 if (nt == 2) {
00666 return;
00667 }
00668 for (i = 3; i <= nt; ++i) {
00669 if (ntot >= MAX_SECONDARIES) {
00670 return;
00671 }
00672 eve[ntot++] = pv[i];
00673 }
00674
00675 }
00676
00677
00678 G4double G4AntiNeutronAnnihilationAtRest::ExNu(G4float ek1)
00679 {
00680 G4float ret_val, r__1;
00681
00682 static G4float cfa, gfa, ran1, ran2, ekin1, atno3;
00683 static G4int magic;
00684 static G4float fpdiv;
00685
00686
00687
00688
00689
00690
00691 ret_val = G4float(0.);
00692 if (targetAtomicMass >= G4float(1.5)) {
00693 magic = 0;
00694 if (G4int(targetCharge + G4float(.1)) == 82) {
00695 magic = 1;
00696 }
00697 ekin1 = ek1;
00698 if (ekin1 < G4float(.1)) {
00699 ekin1 = G4float(.1);
00700 }
00701 if (ekin1 > G4float(4.)) {
00702 ekin1 = G4float(4.);
00703 }
00704
00705
00706 cfa = G4float(.13043478260869565);
00707 cfa = cfa * std::log(ekin1) + G4float(.35);
00708 if (cfa < G4float(.15)) {
00709 cfa = G4float(.15);
00710 }
00711 ret_val = cfa * G4float(7.716) * std::exp(-G4double(cfa));
00712 atno3 = targetAtomicMass;
00713 if (atno3 > G4float(120.)) {
00714 atno3 = G4float(120.);
00715 }
00716 cfa = (atno3 - G4float(1.)) /
00717 G4float(120.) * std::exp(-G4double(atno3 - G4float(1.)) / G4float(120.));
00718 ret_val *= cfa;
00719 r__1 = ekin1;
00720 fpdiv = G4float(1.) - r__1 * r__1 * G4float(.25);
00721 if (fpdiv < G4float(.5)) {
00722 fpdiv = G4float(.5);
00723 }
00724 gfa = (targetAtomicMass - G4float(1.)) /
00725 G4float(70.) * G4float(2.) *
00726 std::exp(-G4double(targetAtomicMass - G4float(1.)) / G4float(70.));
00727 evapEnergy1 = ret_val * fpdiv;
00728 evapEnergy3 = ret_val - evapEnergy1;
00729 Normal(&ran1);
00730 Normal(&ran2);
00731 if (magic == 1) {
00732 ran1 = G4float(0.);
00733 ran2 = G4float(0.);
00734 }
00735 evapEnergy1 *= ran1 * gfa + G4float(1.);
00736 if (evapEnergy1 < G4float(0.)) {
00737 evapEnergy1 = G4float(0.);
00738 }
00739 evapEnergy3 *= ran2 * gfa + G4float(1.);
00740 if (evapEnergy3 < G4float(0.)) {
00741 evapEnergy3 = G4float(0.);
00742 }
00743 while ((ret_val = evapEnergy1 + evapEnergy3) >= ek1) {
00744 evapEnergy1 *= G4float(1.) - G4UniformRand() * G4float(.5);
00745 evapEnergy3 *= G4float(1.) - G4UniformRand() * G4float(.5);
00746 }
00747 }
00748 return ret_val;
00749
00750 }