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