Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4PreCompoundFragment.cc
Go to the documentation of this file.
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 // $Id$
27 //
28 // J. M. Quesada (August 2008).
29 // Based on previous work by V. Lara
30 //
31 // Modified:
32 // 06.09.2008 JMQ Also external choice has been added for:
33 // - superimposed Coulomb barrier (if useSICB=true)
34 // 20.08.2010 V.Ivanchenko cleanup
35 //
36 
37 #include "G4PreCompoundFragment.hh"
38 
39 G4PreCompoundFragment::
40 G4PreCompoundFragment(const G4ParticleDefinition* part,
41  G4VCoulombBarrier* aCoulombBarrier)
42  : G4VPreCompoundFragment(part,aCoulombBarrier)
43 {}
44 
46 {}
47 
50 {
51  //G4cout << theCoulombBarrier << " " << GetMaximalKineticEnergy() << G4endl;
52  // If theCoulombBarrier effect is included in the emission probabilities
53  //if (GetMaximalKineticEnergy() <= 0.0)
54  G4double limit = 0.0;
55  if(OPTxs==0 || useSICB) { limit = theCoulombBarrier; }
56  if (GetMaximalKineticEnergy() <= limit)
57  {
59  return 0.0;
60  }
61  // If theCoulombBarrier effect is included in the emission probabilities
62  // G4double LowerLimit = 0.;
63  // Coulomb barrier is the lower limit
64  // of integration over kinetic energy
65  G4double LowerLimit = limit;
66 
67  // Excitation energy of nucleus after fragment emission is the upper
68  //limit of integration over kinetic energy
69  G4double UpperLimit = GetMaximalKineticEnergy();
70 
72  IntegrateEmissionProbability(LowerLimit,UpperLimit,aFragment);
73  /*
74  G4cout << "## G4PreCompoundFragment::CalcEmisProb "
75  << "Z= " << aFragment.GetZ_asInt()
76  << " A= " << aFragment.GetA_asInt()
77  << " Elow= " << LowerLimit/MeV
78  << " Eup= " << UpperLimit/MeV
79  << " prob= " << theEmissionProbability
80  << G4endl;
81  */
83 }
84 
85 G4double G4PreCompoundFragment::
86 IntegrateEmissionProbability(G4double Low, G4double Up,
87  const G4Fragment & aFragment)
88 {
89  static const G4int N = 10;
90  // 10-Points Gauss-Legendre abcisas and weights
91  static const G4double w[N] = {
92  0.0666713443086881,
93  0.149451349150581,
94  0.219086362515982,
95  0.269266719309996,
96  0.295524224714753,
97  0.295524224714753,
98  0.269266719309996,
99  0.219086362515982,
100  0.149451349150581,
101  0.0666713443086881
102  };
103  static const G4double x[N] = {
104  -0.973906528517172,
105  -0.865063366688985,
106  -0.679409568299024,
107  -0.433395394129247,
108  -0.148874338981631,
109  0.148874338981631,
110  0.433395394129247,
111  0.679409568299024,
112  0.865063366688985,
113  0.973906528517172
114  };
115 
116  G4double Total = 0.0;
117 
118  for (G4int i = 0; i < N; ++i)
119  {
120  G4double KineticE = 0.5*((Up-Low)*x[i]+(Up+Low));
121  Total += w[i]*ProbabilityDistributionFunction(KineticE, aFragment);
122  }
123  Total *= 0.5*(Up-Low);
124  return Total;
125 }
126 
128 GetKineticEnergy(const G4Fragment & aFragment)
129 {
130  //let's keep this way for consistency with CalcEmissionProbability method
131  G4double V = 0.0;
132  if(OPTxs==0 || useSICB) { V = theCoulombBarrier; }
133 
135  if(Tmax < V) { return 0.0; }
136  G4double T(0.0);
137  G4double Probability(1.0);
138  G4double maxProbability = GetEmissionProbability();
139  do
140  {
141  T = V + G4UniformRand()*(Tmax-V);
142  Probability = ProbabilityDistributionFunction(T,aFragment);
143  } while (maxProbability*G4UniformRand() > Probability);
144  return T;
145 }