Geant4  10.01.p03
G4PreCompoundAlpha.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: G4PreCompoundAlpha.cc 90591 2015-06-04 13:45:29Z gcosmo $
27 //
28 // -------------------------------------------------------------------
29 //
30 // GEANT4 Class file
31 //
32 //
33 // File name: G4PreCompoundAlpha
34 //
35 // Author: V.Lara
36 //
37 // Modified:
38 // 21.08.2008 J. M. Quesada add choice of options
39 // 20.08.2010 V.Ivanchenko added G4Pow and G4PreCompoundParameters pointers
40 // use int Z and A and cleanup
41 // 05.07.2013 J.M. Quesada FactorialFactor fixed
42 //
43 
44 #include "G4PreCompoundAlpha.hh"
45 #include "G4SystemOfUnits.hh"
46 #include "G4Alpha.hh"
47 
48 const G4double p0 = 10.95;
49 const G4double p1 = -85.2;
50 const G4double p2 = 1146.;
51 const G4double landa0 = 0.0643;
52 const G4double landa1 = -13.96;
53 const G4double mm0 = 781.2;
54 const G4double mu1 = 0.29;
55 const G4double nu0 = -304.7;
56 const G4double nu1 = -470.0;
57 const G4double nu2 = -8.580;
58 
60  : G4PreCompoundIon(G4Alpha::Alpha(), &theAlphaCoulombBarrier)
61 {
62  ResidualA = GetRestA();
63  ResidualZ = GetRestZ();
64  theA = GetA();
65  theZ = GetZ();
69 }
70 
72 {}
73 
75 {
76  return G4double((N-4)*(P-3)*(N-3)*(P-2))*G4double((N-2)*(P-1)*(N-1)*P)/144.0;
77 }
78 
80 {
81  return 4096.0/G4double(A*A*A);
82 }
83 
85 {
86  G4double rj = 0.0;
87  if(nCharged >=2 && (nParticles-nCharged) >=2 ) {
88  G4double denominator =
89  G4double(nParticles*(nParticles-1)*(nParticles-2)*(nParticles-3));
90  rj = 6.0*nCharged*(nCharged-1)*(nParticles-nCharged)*(nParticles-nCharged-1)
91  /denominator;
92  }
93  return rj;
94 }
95 
97 //J. M. Quesada (Dec 2007-June 2008): New inverse reaction cross sections
98 //OPT=0 Dostrovski's parameterization
99 //OPT=1,2 Chatterjee's paramaterization
100 //OPT=3,4 Kalbach's parameterization
101 //
103 {
104  ResidualA = GetRestA();
105  ResidualZ = GetRestZ();
106  theA = GetA();
107  theZ = GetZ();
111 
112  if (OPTxs==0) { return GetOpt0( K); }
113  else if(OPTxs <= 2) { return GetOpt12( K); }
114  else { return GetOpt34( K); }
115 }
116 
118 {
119  G4double C = 0.0;
120  G4int aZ = theZ + ResidualZ;
121  if (aZ <= 30)
122  {
123  C = 0.10;
124  }
125  else if (aZ <= 50)
126  {
127  C = 0.1 - (aZ-30)*0.001;
128  }
129  else if (aZ < 70)
130  {
131  C = 0.08 - (aZ-50)*0.001;
132  }
133  else
134  {
135  C = 0.06;
136  }
137  return 1.0+C;
138 }
139 
140 //
141 //********************* OPT=1,2 : Chatterjee's cross section ********************
142 //(fitting to cross section from Bechetti & Greenles OM potential)
143 
145 {
146  // JMQ xsec is set constant above limit of validity
147  G4double Kc = std::min(K, 50*MeV);
148  G4double landa ,mu ,nu ,p , Ec,q,r,ji,xs;
149 
150  const G4double delta=1.2;
151 
152  Ec = 1.44*theZ*ResidualZ/(1.5*ResidualAthrd+delta);
153  p = p0 + p1/Ec + p2/(Ec*Ec);
154  landa = landa0*ResidualA + landa1;
155  G4double resmu1 = g4pow->powZ(ResidualA,mu1);
156  mu = mm0*resmu1;
157  nu = resmu1*(nu0 + nu1*Ec + nu2*(Ec*Ec));
158  q = landa - nu/(Ec*Ec) - 2*p*Ec;
159  r = mu + 2*nu/Ec + p*(Ec*Ec);
160 
161  ji=std::max(Kc,Ec);
162  if(Kc < Ec) { xs = p*Kc*Kc + q*Kc + r;}
163  else {xs = p*(Kc - ji)*(Kc - ji) + landa*Kc + mu + nu*(2 - Kc/ji)/ji ;}
164 
165  xs = std::max(xs, 0.0);
166  return xs;
167 }
168 
169 // *********** OPT=3,4 : Kalbach's cross sections (from PRECO code)*************
171 // c ** alpha from huizenga and igo
172 {
173  const G4double flow = 1.e-18;
174  const G4double spill= 1.e+18;
175  const G4double ra = 1.20;
176  const G4double signor = 1.0;
177 
178  //JMQ 13/02/09 increase of reduced radius to lower the barrier
179  // ec = 1.44 * theZ * ResidualZ / (1.5*ResidualAthrd+ra);
180  G4double ec = 1.44 * theZ * ResidualZ / (1.7*ResidualAthrd+ra);
181  G4double ecsq = ec * ec;
182  G4double p = p0 + p1/ec + p2/ecsq;
183  G4double landa = landa0*ResidualA + landa1;
185  G4double mu = mm0 * a;
186  G4double nu = a* (nu0+nu1*ec+nu2*ecsq);
187  G4double xnulam = nu / landa;
188  G4double etest = 0.0;
189  if (xnulam > spill) { xnulam=0.; }
190  else if (xnulam >= flow) { etest = 1.2 *std::sqrt(xnulam); }
191 
192  a = -2.*p*ec + landa - nu/ecsq;
193  G4double b = p*ecsq + mu + 2.*nu/ec;
194  G4double ecut = 0.;
195  G4double cut = a*a - 4.*p*b;
196  if (cut > 0.) { ecut = std::sqrt(cut); }
197  ecut = (ecut-a) / (2*p);
198 
199  //JMQ 290310 for avoiding unphysical increase below minimum (at ecut)
200  // ecut<0 means that there is no cut with energy axis, i.e. xs is set
201  // to 0 bellow minimum
202 
203  G4double elab = K * FragmentA / G4double(ResidualA);
204  G4double sig = 0.;
205 
206  if (elab <= ec) {
207  if (elab > ecut) { sig = std::max(0.0,(p*elab*elab+a*elab+b) * signor); }
208 
209  } else {
210  sig = (landa*elab+mu+nu/elab) * signor;
211  G4double geom = 0.;
212  if (xnulam >= flow && elab >= etest) {
213  geom = std::sqrt(theA*K);
214  geom = 1.23*ResidualAthrd + ra + 4.573/geom;
215  geom = 31.416 * geom * geom;
216  }
217  sig = std::max(sig, geom);
218  }
219  return sig;
220 }
virtual G4double CrossSection(G4double ekin)
static const double MeV
Definition: G4SIunits.hh:193
G4double GetOpt34(G4double K)
G4double ResidualA13() const
G4int GetA() const
const G4double nu1
G4double a
Definition: TRTMaterials.hh:39
G4double GetOpt0(G4double ekin)
int G4int
Definition: G4Types.hh:78
virtual G4double GetAlpha()
const G4double mm0
G4double Z13(G4int Z) const
Definition: G4Pow.hh:127
const G4double landa1
const G4double nu0
const G4double p2
const G4double p1
virtual G4double GetRj(G4int NumberParticles, G4int NumberCharged)
G4int GetRestZ() const
static const G4double A[nN]
const G4double mu1
const G4double p0
T max(const T t1, const T t2)
brief Return the largest of the two arguments
const G4double nu2
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
virtual G4double FactorialFactor(G4int N, G4int P)
G4int GetRestA() const
G4int GetZ() const
virtual G4double CoalescenceFactor(G4int A)
G4double powZ(G4int Z, G4double y) const
Definition: G4Pow.hh:256
double G4double
Definition: G4Types.hh:76
const G4double landa0
G4double GetOpt12(G4double K)