Geant4  10.00.p02
G4NuclearAbrasionGeometry.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 // * *
21 // * Parts of this code which have been developed by QinetiQ Ltd *
22 // * under contract to the European Space Agency (ESA) are the *
23 // * intellectual property of ESA. Rights to use, copy, modify and *
24 // * redistribute this software for general public use are granted *
25 // * in compliance with any licensing, distribution and development *
26 // * policy adopted by the Geant4 Collaboration. This code has been *
27 // * written by QinetiQ Ltd for the European Space Agency, under ESA *
28 // * contract 17191/03/NL/LvH (Aurora Programme). *
29 // * *
30 // * By using, copying, modifying or distributing the software (or *
31 // * any work based on the software) you agree to acknowledge its *
32 // * use in resulting scientific publications, and indicate your *
33 // * acceptance of all terms of the Geant4 Software license. *
34 // ********************************************************************
35 //
36 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
37 //
38 // MODULE: G4NuclearAbrasionGeometry.cc
39 //
40 // Version: B.1
41 // Date: 15/04/04
42 // Author: P R Truscott
43 // Organisation: QinetiQ Ltd, UK
44 // Customer: ESA/ESTEC, NOORDWIJK
45 // Contract: 17191/03/NL/LvH
46 //
47 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
48 //
49 // CHANGE HISTORY
50 // --------------
51 //
52 // 18 November 2003, P R Truscott, QinetiQ Ltd, UK
53 // Created.
54 //
55 // 15 March 2004, P R Truscott, QinetiQ Ltd, UK
56 // Beta release
57 //
58 // 4 June 2004, J.P. Wellisch, CERN, Switzerland
59 // resolving technical portability issues.
60 //
61 // 12 June 2012, A. Ribon, CERN, Switzerland
62 // Fixing trivial warning errors of shadowed variables.
63 //
64 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
66 //
68 #include "G4WilsonRadius.hh"
69 #include "G4PhysicalConstants.hh"
70 #include "G4SystemOfUnits.hh"
72 //
74  G4double AT1, G4double r1)
75 {
76 //
77 //
78 // Initialise variables for interaction geometry.
79 //
80  G4WilsonRadius aR;
81  AP = AP1;
82  AT = AT1;
83  rP = aR.GetWilsonRadius(AP);
84  rT = aR.GetWilsonRadius(AT);
85  r = r1;
86  n = rP / (rP + rT);
87  b = r / (rP + rT);
88  m = rT / rP;
89  Q = (1.0 - b)/n;
90  S = Q * Q;
91  T = S * Q;
92  R = std::sqrt(m*n);
93  U = 1.0/m - 2.0;
94 //
95 //
96 // Initialise the threshold radius-ratio at which interactions are considered
97 // peripheral or central.
98 //
99  rth = 2.0/3.0;
100  B = 10.0 * MeV;
101 }
103 //
105 {;}
107 //
109  {if (rth1 > 0.0 && rth1 <= 1.0) rth = rth1;}
111 //
113  {return rth;}
115 //
117 {
118 //
119 //
120 // Initialise the value for P, then determine the actual value depending upon
121 // whether the projectile is larger or smaller than the target and these radii
122 // in relation to the impact parameter.
123 //
124  G4double valueP = 0.0;
125 
126  if (rT > rP)
127  {
128  if (rT-rP<=r && r<=rT+rP) valueP = 0.125*R*U*S - 0.125*(0.5*R*U+1.0)*T;
129  else valueP = -1.0;
130  }
131  else
132  {
133  if (rP-rT<=r && r<=rP+rT) valueP = 0.125*R*U*S - 0.125*(0.5*std::sqrt(n/m)*U-
134  (std::sqrt(1.0-m*m)/n - 1.0)*std::sqrt((2.0-m)/std::pow(m,5.0)))*T;
135  else valueP = (std::sqrt(1.0-m*m)/n-1.0)*std::sqrt(1.0-b*b/n/n);
136  }
137 
138  if (!(valueP <= 1.0 && valueP>= -1.0))
139  {
140  if (valueP > 1.0) valueP = 1.0;
141  else valueP = -1.0;
142  }
143  return valueP;
144 }
146 //
148 {
149 //
150 //
151 // Initialise the value for F, then determine the actual value depending upon
152 // whether the projectile is larger or smaller than the target and these radii
153 // in relation to the impact parameter.
154 //
155  G4double valueF = 0.0;
156 
157  if (rT > rP)
158  {
159  if (rT-rP<=r && r<=rT+rP) valueF = 0.75*R*S - 0.125*(3.0*R-1.0)*T;
160  else valueF = 1.0;
161  }
162  else
163  {
164  if (rP-rT<=r && r<=rP+rT) valueF = 0.75*R*S - 0.125*(3.0*std::sqrt(n/m)-
165  (1.0-std::pow(1.0-m*m,3.0/2.0))*std::sqrt(1.0-std::pow(1.0-m,2.0))/std::pow(m,3.0))*T;
166  else valueF = (1.0-std::pow(1.0-m*m,3.0/2.0))*std::sqrt(1.0-b*b/n/n);
167  }
168 
169  if (!(valueF <= 1.0 && valueF>= 0.0))
170  {
171  if (valueF > 1.0) valueF = 1.0;
172  else valueF = 0.0;
173  }
174  return valueF;
175 }
177 //
179 {
180  G4double F1 = F();
181  G4double P1 = P();
182  G4double Es = 0.0;
183 
184  Es = 0.95 * MeV * 4.0 * pi * rP*rP/fermi/fermi *
185  (1.0+P1-std::pow(1.0-F1,2.0/3.0));
186 // if (rT < rP && r < rP-rT)
187  if ((r-rP)/rT < rth)
188  {
189  G4double omega = 0.0;
190  if (AP < 12.0) omega = 1500.0;
191  else if (AP <= 16.0) omega = 1500.0 - 320.0*(AP-12.0);
192  Es *= 1.0 + F1*(5.0+omega*F1*F1);
193  }
194 
195  if (Es < 0.0)
196  Es = 0.0;
197  else if (Es > B * AP)
198  Es = B * AP;
199  return Es;
200 }
201 
202 
204 {
205  // This member function declares a new G4NuclearAbrasionGeometry object
206  // but with the projectile and target exchanged to determine the values
207  // for F and P. Determination of the excess surface area and excitation
208  // energy is as above.
209 
210  G4NuclearAbrasionGeometry* revAbrasionGeometry =
212  G4double F1 = revAbrasionGeometry->F();
213  G4double P1 = revAbrasionGeometry->P();
214  G4double Es = 0.0;
215 
216  Es = 0.95 * MeV * 4.0 * pi * rT*rT/fermi/fermi *
217  (1.0+P1-std::pow(1.0-F1,2.0/3.0));
218 
219 // if (rP < rT && r < rT-rP)
220  if ((r-rT)/rP < rth) {
221  G4double omega = 0.0;
222  if (AT < 12.0) omega = 1500.0;
223  else if (AT <= 16.0) omega = 1500.0 - 320.0*(AT-12.0);
224  Es *= 1.0 + F1*(5.0+omega*F1*F1);
225  }
226 
227  if (Es < 0.0)
228  Es = 0.0;
229  else if (Es > B * AT)
230  Es = B * AT;
231 
232  delete revAbrasionGeometry;
233 
234  return Es;
235 }
static const double MeV
Definition: G4SIunits.hh:193
G4NuclearAbrasionGeometry(G4double AP, G4double AT, G4double r)
static const G4double * P1[nN]
const G4double pi
double G4double
Definition: G4Types.hh:76
static const double fermi
Definition: G4SIunits.hh:93
G4double GetWilsonRadius(G4double A)