source: trunk/source/processes/hadronic/stopping/src/G4AntiNeutronAnnihilationAtRest.cc @ 1337

Last change on this file since 1337 was 1055, checked in by garnier, 15 years ago

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