Geant4  10.02.p01
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: G4PreCompoundFragment.cc 94281 2015-11-10 15:00:15Z gcosmo $
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 #include "G4KalbachCrossSection.hh"
40 #include "Randomize.hh"
41 
42 const G4int NPOINTSGL = 10;
43 // 10-Points Gauss-Legendre abcisas and weights
44 const G4double w[NPOINTSGL] = {
45  0.0666713443086881,
46  0.149451349150581,
47  0.219086362515982,
48  0.269266719309996,
49  0.295524224714753,
50  0.295524224714753,
51  0.269266719309996,
52  0.219086362515982,
53  0.149451349150581,
54  0.0666713443086881
55  };
56 const G4double x[NPOINTSGL] = {
57  -0.973906528517172,
58  -0.865063366688985,
59  -0.679409568299024,
60  -0.433395394129247,
61  -0.148874338981631,
62  0.148874338981631,
63  0.433395394129247,
64  0.679409568299024,
65  0.865063366688985,
66  0.973906528517172
67 };
68 
71  G4VCoulombBarrier* aCoulombBarrier)
72  : G4VPreCompoundFragment(part,aCoulombBarrier)
73 {
74  muu = probmax = 0.0;
75  index = 0;
76  if(1 == theZ) { index = theA; }
77  else { index = theA + 1; }
78  for(G4int i=0; i<12; ++i) { probability[i] = 0.0; }
79 }
80 
82 {}
83 
86 {
87  //G4cout << theCoulombBarrier << " " << GetMaximalKineticEnergy() << G4endl;
88  // If theCoulombBarrier effect is included in the emission probabilities
89  // Coulomb barrier is the lower limit of integration over kinetic energy
90 
91  G4double LowerLimit = 0.0;
92  if(OPTxs==0 || useSICB) { LowerLimit = theCoulombBarrier; }
93 
94  if (theMaxKinEnergy <= LowerLimit)
95  {
97  return 0.0;
98  }
99 
100  // compute power once
101  if(OPTxs <= 2) {
103  } else {
105  }
106 
108  IntegrateEmissionProbability(LowerLimit,theMaxKinEnergy,aFragment);
109  /*
110  G4cout << "## G4PreCompoundFragment::CalcEmisProb "
111  << "Z= " << aFragment.GetZ_asInt()
112  << " A= " << aFragment.GetA_asInt()
113  << " Elow= " << LowerLimit/MeV
114  << " Eup= " << UpperLimit/MeV
115  << " prob= " << theEmissionProbability
116  << G4endl;
117  */
118  return theEmissionProbability;
119 }
120 
123  const G4Fragment & aFragment)
124 {
125  G4double sum = 0.0;
126  G4double del = (up - low)*0.5;
127  G4double avr = (up + low)*0.5;
128 
129  G4double e, y;
130  probmax = 0.0;
131 
132  for (G4int i=0; i<NPOINTSGL; ++i) {
133  e = del*x[i] + avr;
134  y = ProbabilityDistributionFunction(e, aFragment);
135  probability[i] = y;
136  probmax = std::max(probmax, y);
137  sum += w[i]*y;
138  }
139  return sum*del;
140 }
141 
143 {
144  G4double res;
145  if(OPTxs == 0) {
146  res = GetOpt0(ekin);
147  } else if(OPTxs <= 2) {
149  index, theZ,
150  theResZ, theResA);
151  } else {
153  index, theZ, theA,
154  theResZ, theResA);
155  }
156  return res;
157 }
158 
159 // *********************** OPT=0 : Dostrovski's cross section ***************
161 {
163  // cross section is now given in mb (r0 is in mm) for the sake of consistency
164  //with the rest of the options
165  return 1.e+25*CLHEP::pi*r0*r0*theResA13*GetAlpha()*(1.+GetBeta()/ekin);
166 }
167 
169 {
170  //let's keep this way for consistency with CalcEmissionProbability method
171  G4double limit = 0.0;
172  if(useSICB) { limit = theCoulombBarrier; }
173 
174  if(theMaxKinEnergy <= limit) { return 0.0; }
175 
176  G4double delta = theMaxKinEnergy - limit;
177  static const G4double toler = 1.25;
178  probmax *= toler;
179  G4double prob, T(0.0);
180  CLHEP::HepRandomEngine* rndm = G4Random::getTheEngine();
181  G4int i;
182  for(i=0; i<100; ++i) {
183  T = limit + delta*rndm->flat();
184  prob = ProbabilityDistributionFunction(T, fragment);
185  /*
186  if(prob > probmax) {
187  G4cout << "G4PreCompoundFragment WARNING: prob= " << prob
188  << " probmax= " << probmax << G4endl;
189  G4cout << "i= " << i << " Z= " << theZ << " A= " << theA
190  << " resZ= " << theResZ << " resA= " << theResA << "\n"
191  << " T= " << T << " Tmax= " << theMaxKinEnergy
192  << " Tmin= " << limit
193  << G4endl;
194  for(G4int i=0; i<N; ++i) { G4cout << " " << probability[i]; }
195  G4cout << G4endl;
196  }
197  */
198  // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
199  if(probmax*rndm->flat() <= prob) { break; }
200  }
201  //G4cout << "G4PreCompoundFragment: i= " << i << " T(MeV)= " << T
202  //<< " Emin(MeV)= " << limit << " Emax= " << theMaxKinEnergy << G4endl;
203  return T;
204 }
205 
G4double GetOpt0(G4double ekin) const
G4double SampleKineticEnergy(const G4Fragment &aFragment)
const G4double w[NPOINTSGL]
static G4double ComputePowerParameter(G4int resA, G4int idx)
static G4double ComputePowerParameter(G4int resA, G4int idx)
static G4double ComputeCrossSection(G4double K, G4double resA13, G4double amu1, G4int idx, G4int Z, G4int resZ, G4int resA)
int G4int
Definition: G4Types.hh:78
virtual G4double GetAlpha() const =0
virtual G4double GetBeta() const =0
G4double CrossSection(G4double ekin) const
G4double CalcEmissionProbability(const G4Fragment &aFragment)
const G4int NPOINTSGL
static const double pi
Definition: G4SIunits.hh:74
T max(const T t1, const T t2)
brief Return the largest of the two arguments
virtual G4double ProbabilityDistributionFunction(G4double K, const G4Fragment &aFragment)=0
const G4double x[NPOINTSGL]
G4PreCompoundParameters * theParameters
G4double IntegrateEmissionProbability(G4double Low, G4double Up, const G4Fragment &aFragment)
double G4double
Definition: G4Types.hh:76
static G4double ComputeCrossSection(G4double K, G4double resA13, G4double amu1, G4int idx, G4int Z, G4int A, G4int resZ, G4int resA)
const G4double r0