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

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

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