ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4AntiNeutronAnnihilationAtRest.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4AntiNeutronAnnihilationAtRest.cc
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // G4AntiNeutronAnnihilationAtRest physics process
27 // Larry Felawka (TRIUMF), April 1998
28 //---------------------------------------------------------------------
29 
31 #include "G4SystemOfUnits.hh"
32 #include "G4DynamicParticle.hh"
33 #include "G4ParticleTypes.hh"
35 #include "G4HadronicDeprecate.hh"
36 #include "Randomize.hh"
37 #include "G4Exp.hh"
38 #include "G4Log.hh"
39 #include "G4Pow.hh"
40 
41 #define MAX_SECONDARIES 100
42 
43 // constructor
44 
46  G4ProcessType aType) :
47  G4VRestProcess (processName, aType), // initialization
48  massPionMinus(G4PionMinus::PionMinus()->GetPDGMass()/GeV),
49  massPionZero(G4PionZero::PionZero()->GetPDGMass()/GeV),
50  massPionPlus(G4PionPlus::PionPlus()->GetPDGMass()/GeV),
51  massGamma(G4Gamma::Gamma()->GetPDGMass()/GeV),
52  massAntiNeutron(G4AntiNeutron::AntiNeutron()->GetPDGMass()/GeV),
53  massNeutron(G4Neutron::Neutron()->GetPDGMass()/GeV),
54  pdefGamma(G4Gamma::Gamma()),
55  pdefPionPlus(G4PionPlus::PionPlus()),
56  pdefPionZero(G4PionZero::PionZero()),
57  pdefPionMinus(G4PionMinus::PionMinus()),
58  pdefProton(G4Proton::Proton()),
59  pdefNeutron(G4Neutron::Neutron()),
60  pdefAntiNeutron(G4AntiNeutron::AntiNeutron()),
61  pdefDeuteron(G4Deuteron::Deuteron()),
62  pdefTriton(G4Triton::Triton()),
63  pdefAlpha(G4Alpha::Alpha())
64 {
65  G4HadronicDeprecate("G4AntiNeutronAnnihilationAtRest");
66  if (verboseLevel>0) {
67  G4cout << GetProcessName() << " is created "<< G4endl;
68  }
73 
76  = evapEnergy3 = 0.0;
77  ngkine = ntot = 0;
78 }
79 
80 // destructor
81 
83 {
85  delete [] pv;
86  delete [] eve;
87  delete [] gkin;
88 }
89 
91 {
93 }
94 
96 {
98 }
99 
100 // methods.............................................................................
101 
104  )
105 {
106  return ( &particle == pdefAntiNeutron );
107 
108 }
109 
110 // Warning - this method may be optimized away if made "inline"
112 {
113  return ( ngkine );
114 
115 }
116 
117 // Warning - this method may be optimized away if made "inline"
119 {
120  return ( &gkin[0] );
121 
122 }
123 
125  const G4Track& track,
127  )
128 {
129  // beggining of tracking
131 
132  // condition is set to "Not Forced"
133  *condition = NotForced;
134 
135  // get mean life time
136  currentInteractionLength = GetMeanLifeTime(track, condition);
137 
138  if ((currentInteractionLength <0.0) || (verboseLevel>2)){
139  G4cout << "G4AntiNeutronAnnihilationAtRestProcess::AtRestGetPhysicalInteractionLength ";
140  G4cout << "[ " << GetProcessName() << "]" <<G4endl;
141  track.GetDynamicParticle()->DumpInfo();
142  G4cout << " in Material " << track.GetMaterial()->GetName() <<G4endl;
143  G4cout << "MeanLifeTime = " << currentInteractionLength/ns << "[ns]" <<G4endl;
144  }
145 
147 
148 }
149 
151  const G4Track& track,
152  const G4Step&
153  )
154 //
155 // Handles AntiNeutrons at rest; an AntiNeutron can either create secondaries
156 // or do nothing (in which case it should be sent back to decay-handling
157 // section
158 //
159 {
160 
161 // Initialize ParticleChange
162 // all members of G4VParticleChange are set to equal to
163 // corresponding member in G4Track
164 
166 
167 // Store some global quantities that depend on current material and particle
168 
169  globalTime = track.GetGlobalTime()/s;
170  G4Material * aMaterial = track.GetMaterial();
171  const G4int numberOfElements = aMaterial->GetNumberOfElements();
172  const G4ElementVector* theElementVector = aMaterial->GetElementVector();
173 
174  const G4double* theAtomicNumberDensity = aMaterial->GetAtomicNumDensityVector();
175  G4double normalization = 0;
176  for ( G4int i1=0; i1 < numberOfElements; i1++ )
177  {
178  normalization += theAtomicNumberDensity[i1] ; // change when nucleon specific
179  // probabilities are included.
180  }
181  G4double runningSum= 0.;
182  G4double random = G4UniformRand()*normalization;
183  for ( G4int i2=0; i2 < numberOfElements; i2++ )
184  {
185  runningSum += theAtomicNumberDensity[i2]; // change when nucleon specific
186  // probabilities are included.
187  if (random<=runningSum)
188  {
189  targetCharge = G4double( ((*theElementVector)[i2])->GetZ());
190  targetAtomicMass = (*theElementVector)[i2]->GetN();
191  }
192  }
193  if (random>runningSum)
194  {
195  targetCharge = G4double( ((*theElementVector)[numberOfElements-1])->GetZ());
196  targetAtomicMass = (*theElementVector)[numberOfElements-1]->GetN();
197  }
198 
199  if (verboseLevel>1) {
200  G4cout << "G4AntiNeutronAnnihilationAtRest::AtRestDoIt is invoked " <<G4endl;
201  }
202 
204  G4float localtime;
205 
207 
208  GenerateSecondaries(); // Generate secondaries
209 
211 
212  for ( G4int isec = 0; isec < ngkine; isec++ ) {
213  G4DynamicParticle* aNewParticle = new G4DynamicParticle;
214  aNewParticle->SetDefinition( gkin[isec].GetParticleDef() );
215  aNewParticle->SetMomentum( gkin[isec].GetMomentum() * GeV );
216 
217  localtime = globalTime + gkin[isec].GetTOF();
218 
219  G4Track* aNewTrack = new G4Track( aNewParticle, localtime*s, position );
220  aNewTrack->SetTouchableHandle(track.GetTouchableHandle());
221  aParticleChange.AddSecondary( aNewTrack );
222 
223  }
224 
226 
227  aParticleChange.ProposeTrackStatus(fStopAndKill); // Kill the incident AntiNeutron
228 
229 // clear InteractionLengthLeft
230 
232 
233  return &aParticleChange;
234 
235 }
236 
237 
239 {
240  G4int index;
241  G4int l;
242  G4int nopt;
243  G4int i;
244  // DHW 15 May 2011: unused: static G4ParticleDefinition* jnd;
245 
246  for (i = 1; i <= MAX_SECONDARIES; ++i) {
247  pv[i].SetZero();
248  }
249 
250 
251  ngkine = 0; // number of generated secondary particles
252  ntot = 0;
253  result.SetZero();
256  result.SetTOF( 0. );
258 
259  // *** SELECT PROCESS FOR CURRENT PARTICLE ***
260 
262 
263  // *** CHECK WHETHER THERE ARE NEW PARTICLES GENERATED ***
264  if (ntot != 0 || result.GetParticleDef() != pdefAntiNeutron) {
265  // *** CURRENT PARTICLE IS NOT THE SAME AS IN THE BEGINNING OR/AND ***
266  // *** ONE OR MORE SECONDARIES HAVE BEEN GENERATED ***
267 
268  // --- INITIAL PARTICLE TYPE HAS BEEN CHANGED ==> PUT NEW TYPE ON ---
269  // --- THE GEANT TEMPORARY STACK ---
270 
271  // --- PUT PARTICLE ON THE STACK ---
272  gkin[0] = result;
273  gkin[0].SetTOF( result.GetTOF() * 5e-11 );
274  ngkine = 1;
275 
276  // --- ALL QUANTITIES ARE TAKEN FROM THE GHEISHA STACK WHERE THE ---
277  // --- CONVENTION IS THE FOLLOWING ---
278 
279  // --- ONE OR MORE SECONDARIES HAVE BEEN GENERATED ---
280  for (l = 1; l <= ntot; ++l) {
281  index = l - 1;
282  // DHW 15 May 2011: unused: jnd = eve[index].GetParticleDef();
283 
284  // --- ADD PARTICLE TO THE STACK IF STACK NOT YET FULL ---
285  if (ngkine < MAX_SECONDARIES) {
286  gkin[ngkine] = eve[index];
287  gkin[ngkine].SetTOF( eve[index].GetTOF() * 5e-11 );
288  ++ngkine;
289  }
290  }
291  }
292  else {
293  // --- NO SECONDARIES GENERATED AND PARTICLE IS STILL THE SAME ---
294  // --- ==> COPY EVERYTHING BACK IN THE CURRENT GEANT STACK ---
295  ngkine = 0;
296  ntot = 0;
297  globalTime += result.GetTOF() * G4float(5e-11);
298  }
299 
300  // --- LIMIT THE VALUE OF NGKINE IN CASE OF OVERFLOW ---
302 
303 } // GenerateSecondaries
304 
305 
307 {
308  G4int i;
309  G4float r, p1, p2, p3;
310  G4int fivex;
311  G4float rr, ran, rrr, ran1;
312 
313  // *** GENERATION OF POISSON DISTRIBUTION ***
314  // *** NVE 16-MAR-1988 CERN GENEVA ***
315  // ORIGIN : H.FESEFELDT (27-OCT-1983)
316 
317  // --- USE NORMAL DISTRIBUTION FOR <X> > 9.9 ---
318  if (xav > G4float(9.9)) {
319  // ** NORMAL DISTRIBUTION WITH SIGMA**2 = <X>
320  Normal(&ran1);
321  ran1 = xav + ran1 * std::sqrt(xav);
322  *iran = G4int(ran1);
323  if (*iran < 0) {
324  *iran = 0;
325  }
326  }
327  else {
328  fivex = G4int(xav * G4float(5.));
329  *iran = 0;
330  if (fivex > 0) {
331  r = G4Exp(-G4double(xav));
332  ran1 = G4UniformRand();
333  if (ran1 > r) {
334  rr = r;
335  for (i = 1; i <= fivex; ++i) {
336  ++(*iran);
337  if (i <= 5) {
338  rrr = G4Pow::GetInstance()->powN(xav, i) / NFac(i);
339  }
340  // ** STIRLING' S FORMULA FOR LARGE NUMBERS
341  if (i > 5) {
342  rrr = G4Exp(i * G4Log(xav) -
343  (i + G4float(.5)) * G4Log(i * G4float(1.)) +
344  i - G4float(.9189385));
345  }
346  rr += r * rrr;
347  if (ran1 <= rr) {
348  break;
349  }
350  }
351  }
352  }
353  else {
354  // ** FOR VERY SMALL XAV TRY IRAN=1,2,3
355  p1 = xav * G4Exp(-G4double(xav));
356  p2 = xav * p1 / G4float(2.);
357  p3 = xav * p2 / G4float(3.);
358  ran = G4UniformRand();
359  if (ran >= p3) {
360  if (ran >= p2) {
361  if (ran >= p1) {
362  *iran = 0;
363  }
364  else {
365  *iran = 1;
366  }
367  }
368  else {
369  *iran = 2;
370  }
371  }
372  else {
373  *iran = 3;
374  }
375  }
376  }
377 
378 } // Poisso
379 
380 
382 {
383  G4int ret_val;
384 
385  G4int i, j;
386 
387  // *** NVE 16-MAR-1988 CERN GENEVA ***
388  // ORIGIN : H.FESEFELDT (27-OCT-1983)
389 
390  ret_val = 1;
391  j = n;
392  if (j > 1) {
393  if (j > 10) {
394  j = 10;
395  }
396  for (i = 2; i <= j; ++i) {
397  ret_val *= i;
398  }
399  }
400  return ret_val;
401 
402 } // NFac
403 
404 
406 {
407  // *** NVE 14-APR-1988 CERN GENEVA ***
408  // ORIGIN : H.FESEFELDT (27-OCT-1983)
409 
410  *ran = (G4float)(-6. + 12.*G4UniformRand());
411 } // Normal
412 
413 
415 {
416  G4float brr[3] = { G4float(.125),G4float(.25),G4float(.5) };
417 
418  G4float r__1;
419 
420  G4int i, ii, kk;
421  G4int nt;
422  G4float cfa, eka;
423  G4int ika, nbl;
424  G4float ran, pcm;
425  G4int isw;
426  G4float tex;
427  G4ParticleDefinition* ipa1;
428  G4float ran1, ran2, ekin, tkin;
429  G4float targ;
431  G4float ekin1, ekin2, black;
432  G4float pnrat, rmnve1, rmnve2;
433  G4float ek, en;
434 
435  // *** ANTI NEUTRON ANNIHILATION AT REST ***
436  // *** NVE 04-MAR-1988 CERN GENEVA ***
437  // ORIGIN : H.FESEFELDT (09-JULY-1987)
438 
439  // NOPT=0 NO ANNIHILATION
440  // NOPT=1 ANNIH.IN PI+ PI-
441  // NOPT=2 ANNIH.IN PI0 PI0
442  // NOPT=3 ANNIH.IN PI+ PI0
443  // NOPT=4 ANNIH.IN GAMMA GAMMA
444 
445  pv[1].SetZero();
446  pv[1].SetMass( massAntiNeutron );
447  pv[1].SetKineticEnergyAndUpdate( 0. );
448  pv[1].SetTOF( result.GetTOF() );
450  isw = 1;
451  ran = G4UniformRand();
452  if (ran > brr[0]) {
453  isw = 2;
454  }
455  if (ran > brr[1]) {
456  isw = 3;
457  }
458  if (ran > brr[2]) {
459  isw = 4;
460  }
461  *nopt = isw;
462  // **
463  // ** EVAPORATION
464  // **
465  rmnve1 = massPionPlus;
466  rmnve2 = massPionMinus;
467  if (isw == 2) {
468  rmnve1 = massPionZero;
469  rmnve2 = massPionZero;
470  }
471  if (isw == 3) {
472  rmnve2 = massPionZero;
473  }
474  if (isw == 4) {
475  rmnve1 = massGamma;
476  rmnve2 = massGamma;
477  }
478  ek = massNeutron + massAntiNeutron - rmnve1 - rmnve2;
479  tkin = ExNu(ek);
480  ek -= tkin;
481  if (ek < G4float(1e-4)) {
482  ek = G4float(1e-4);
483  }
484  ek /= G4float(2.);
485  en = ek + (rmnve1 + rmnve2) / G4float(2.);
486  r__1 = en * en - rmnve1 * rmnve2;
487  pcm = r__1 > 0 ? std::sqrt(r__1) : 0;
488  pv[2].SetZero();
489  pv[2].SetMass( rmnve1 );
490  pv[3].SetZero();
491  pv[3].SetMass( rmnve2 );
492  if (isw > 3) {
493  pv[2].SetMass( 0. );
494  pv[3].SetMass( 0. );
495  }
496  pv[2].SetEnergyAndUpdate( std::sqrt(pv[2].GetMass()*pv[2].GetMass()+pcm*pcm) );
497  pv[2].SetTOF( result.GetTOF() );
498  pv[3].SetEnergy( std::sqrt(pv[3].GetMass()*pv[3].GetMass()+pcm*pcm) );
499  pv[3].SetMomentumAndUpdate( -pv[2].GetMomentum().x(), -pv[2].GetMomentum().y(), -pv[2].GetMomentum().z() );
500  pv[3].SetTOF( result.GetTOF() );
501  switch ((int)isw) {
502  case 1:
505  break;
506  case 2:
509  break;
510  case 3:
513  break;
514  case 4:
517  break;
518  }
519  nt = 3;
520  if (targetAtomicMass >= G4float(1.5)) {
521  cfa = (targetAtomicMass - G4float(1.)) / G4float(120.) *
522  G4float(.025) * G4Exp(-G4double(targetAtomicMass - G4float(1.)) /
523  G4float(120.));
524  targ = G4float(1.);
525  tex = evapEnergy1;
526  if (tex >= G4float(.001)) {
527  black = (targ * G4float(1.25) +
529  Poisso(black, &nbl);
530  if (G4float(G4int(targ) + nbl) > targetAtomicMass) {
531  nbl = G4int(targetAtomicMass - targ);
532  }
533  if (nt + nbl > (MAX_SECONDARIES - 2)) {
534  nbl = (MAX_SECONDARIES - 2) - nt;
535  }
536  if (nbl > 0) {
537  ekin = tex / nbl;
538  ekin2 = 0.0f;
539  for (i = 1; i <= nbl; ++i) {
540  if (nt == (MAX_SECONDARIES - 2)) {
541  continue;
542  }
543  if (ekin2 > tex) {
544  break;
545  }
546  ran1 = G4UniformRand();
547  Normal(&ran2);
548  ekin1 = -G4double(ekin) * G4Log(ran1) -
549  cfa * (ran2 * G4float(.5) + G4float(1.));
550  if (ekin1 < 0.0f) {
551  ekin1 = G4Log(ran1) * G4float(-.01);
552  }
553  ekin1 *= G4float(1.);
554  ekin2 += ekin1;
555  if (ekin2 > tex) {
556  ekin1 = tex - (ekin2 - ekin1);
557  }
558  if (ekin1 < 0.0f) {
559  ekin1 = G4float(.001);
560  }
561  ipa1 = pdefNeutron;
562  pnrat = G4float(1.) - targetCharge / targetAtomicMass;
563  if (G4UniformRand() > pnrat) {
564  ipa1 = pdefProton;
565  }
566  ++nt;
567  pv[nt].SetZero();
568  pv[nt].SetMass( ipa1->GetPDGMass()/GeV );
569  pv[nt].SetKineticEnergyAndUpdate( ekin1 );
570  pv[nt].SetTOF( result.GetTOF() );
571  pv[nt].SetParticleDef( ipa1 );
572  }
573  if (targetAtomicMass >= G4float(230.) && ek <= G4float(2.)) {
574  ii = nt + 1;
575  kk = 0;
576  eka = ek;
577  if (eka > G4float(1.)) {
578  eka *= eka;
579  }
580  if (eka < 0.1f) {
581  eka = 0.1f;
582  }
583  ika = G4int(G4float(3.6) / eka);
584  for (i = 1; i <= nt; ++i) {
585  --ii;
586  if (pv[ii].GetParticleDef() != pdefProton) {
587  continue;
588  }
589  ipa1 = pdefNeutron;
590  pv[ii].SetMass( ipa1->GetPDGMass()/GeV );
591  pv[ii].SetParticleDef( ipa1 );
592  ++kk;
593  if (kk > ika) {
594  break;
595  }
596  }
597  }
598  }
599  }
600  // **
601  // ** THEN ALSO DEUTERONS, TRITONS AND ALPHAS
602  // **
603  tex = evapEnergy3;
604  if (tex >= G4float(.001)) {
605  black = (targ * G4float(1.25) + G4float(1.5)) * evapEnergy3 /
607  Poisso(black, &nbl);
608  if (nt + nbl > (MAX_SECONDARIES - 2)) {
609  nbl = (MAX_SECONDARIES - 2) - nt;
610  }
611  if (nbl > 0) {
612  ekin = tex / nbl;
613  ekin2 = 0.0f;
614  for (i = 1; i <= nbl; ++i) {
615  if (nt == (MAX_SECONDARIES - 2)) {
616  continue;
617  }
618  if (ekin2 > tex) {
619  break;
620  }
621  ran1 = G4UniformRand();
622  Normal(&ran2);
623  ekin1 = -G4double(ekin) * G4Log(ran1) -
624  cfa * (ran2 * G4float(.5) + G4float(1.));
625  if (ekin1 < 0.0f) {
626  ekin1 = G4Log(ran1) * G4float(-.01);
627  }
628  ekin1 *= G4float(1.);
629  ekin2 += ekin1;
630  if (ekin2 > tex) {
631  ekin1 = tex - (ekin2 - ekin1);
632  }
633  if (ekin1 < 0.0f) {
634  ekin1 = G4float(.001);
635  }
636  ran = G4UniformRand();
637  inve = pdefDeuteron;
638  if (ran > G4float(.6)) {
639  inve = pdefTriton;
640  }
641  if (ran > G4float(.9)) {
642  inve = pdefAlpha;
643  }
644  ++nt;
645  pv[nt].SetZero();
646  pv[nt].SetMass( inve->GetPDGMass()/GeV );
647  pv[nt].SetKineticEnergyAndUpdate( ekin1 );
648  pv[nt].SetTOF( result.GetTOF() );
649  pv[nt].SetParticleDef( inve );
650  }
651  }
652  }
653  }
654  result = pv[2];
655  if (nt == 2) {
656  return;
657  }
658  for (i = 3; i <= nt; ++i) {
659  if (ntot >= MAX_SECONDARIES) {
660  return;
661  }
662  eve[ntot++] = pv[i];
663  }
664 
665 } // AntiNeutronAnnihilation
666 
667 
669 {
670  G4float ret_val, r__1;
671 
672  G4float cfa, gfa, ran1, ran2, ekin1, atno3;
673  G4int magic;
674  G4float fpdiv;
675 
676  // *** NUCLEAR EVAPORATION AS FUNCTION OF ATOMIC NUMBER ATNO ***
677  // *** AND KINETIC ENERGY EKIN OF PRIMARY PARTICLE ***
678  // *** NVE 04-MAR-1988 CERN GENEVA ***
679  // ORIGIN : H.FESEFELDT (10-DEC-1986)
680 
681  ret_val = 0.f;
682  if (targetAtomicMass >= G4float(1.5)) {
683  magic = 0;
684  if (G4int(targetCharge + 0.1f) == 82) {
685  magic = 1;
686  }
687  ekin1 = ek1;
688  if (ekin1 < 0.1f) {
689  ekin1 = 0.1f;
690  }
691  if (ekin1 > G4float(4.)) {
692  ekin1 = G4float(4.);
693  }
694  // ** 0.35 VALUE AT 1 GEV
695  // ** 0.05 VALUE AT 0.1 GEV
696  cfa = G4float(.13043478260869565);
697  cfa = cfa * G4Log(ekin1) + G4float(.35);
698  if (cfa < G4float(.15)) {
699  cfa = G4float(.15);
700  }
701  ret_val = cfa * G4float(7.716) * G4Exp(-G4double(cfa));
702  atno3 = targetAtomicMass;
703  if (atno3 > G4float(120.)) {
704  atno3 = G4float(120.);
705  }
706  cfa = (atno3 - G4float(1.)) /
707  G4float(120.) * G4Exp(-G4double(atno3 - G4float(1.)) / G4float(120.));
708  ret_val *= cfa;
709  r__1 = ekin1;
710  fpdiv = G4float(1.) - r__1 * r__1 * G4float(.25);
711  if (fpdiv < G4float(.5)) {
712  fpdiv = G4float(.5);
713  }
714  gfa = (targetAtomicMass - G4float(1.)) /
715  G4float(70.) * G4float(2.) *
717  evapEnergy1 = ret_val * fpdiv;
718  evapEnergy3 = ret_val - evapEnergy1;
719  Normal(&ran1);
720  Normal(&ran2);
721  if (magic == 1) {
722  ran1 = 0.0f;
723  ran2 = 0.0f;
724  }
725  evapEnergy1 *= ran1 * gfa + G4float(1.);
726  if (evapEnergy1 < 0.0f) {
727  evapEnergy1 = 0.0f;
728  }
729  evapEnergy3 *= ran2 * gfa + G4float(1.);
730  if (evapEnergy3 < 0.0f) {
731  evapEnergy3 = 0.0f;
732  }
733 
734  // Loop checking, 06-Aug-2015, Vladimir Ivanchenko
735  while ((ret_val = evapEnergy1 + evapEnergy3) >= ek1) {
736  evapEnergy1 *= G4float(1.) - G4UniformRand() * G4float(.5);
737  evapEnergy3 *= G4float(1.) - G4UniformRand() * G4float(.5);
738  }
739  }
740  return ret_val;
741 
742 } // ExNu