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