source: trunk/source/processes/hadronic/stopping/src/G4AntiProtonAnnihilationAtRest.cc @ 846

Last change on this file since 846 was 819, checked in by garnier, 16 years ago

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