source: trunk/source/processes/hadronic/models/abrasion/src/G4NuclearAbrasionGeometry.cc@ 1036

Last change on this file since 1036 was 819, checked in by garnier, 17 years ago

import all except CVS
File size: 7.7 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// * *
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// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
62////////////////////////////////////////////////////////////////////////////////
63//
64#include "G4NuclearAbrasionGeometry.hh"
65#include "G4WilsonRadius.hh"
66////////////////////////////////////////////////////////////////////////////////
67//
68G4NuclearAbrasionGeometry::G4NuclearAbrasionGeometry (G4double AP1,
69 G4double AT1, G4double r1)
70{
71//
72//
73// Initialise variables for interaction geometry.
74//
75 G4WilsonRadius aR;
76 AP = AP1;
77 AT = AT1;
78 rP = aR.GetWilsonRadius(AP);
79 rT = aR.GetWilsonRadius(AT);
80 r = r1;
81 n = rP / (rP + rT);
82 b = r / (rP + rT);
83 m = rT / rP;
84 Q = (1.0 - b)/n;
85 S = Q * Q;
86 T = S * Q;
87 R = std::sqrt(m*n);
88 U = 1.0/m - 2.0;
89//
90//
91// Initialise the threshold radius-ratio at which interactions are considered
92// peripheral or central.
93//
94 rth = 2.0/3.0;
95 B = 10.0 * MeV;
96}
97////////////////////////////////////////////////////////////////////////////////
98//
99G4NuclearAbrasionGeometry::~G4NuclearAbrasionGeometry ()
100{;}
101////////////////////////////////////////////////////////////////////////////////
102//
103void G4NuclearAbrasionGeometry::SetPeripheralThreshold (G4double rth1)
104 {if (rth1 > 0.0 && rth1 <= 1.0) rth = rth1;}
105////////////////////////////////////////////////////////////////////////////////
106//
107G4double G4NuclearAbrasionGeometry::GetPeripheralThreshold ()
108 {return rth;}
109////////////////////////////////////////////////////////////////////////////////
110//
111G4double G4NuclearAbrasionGeometry::P ()
112{
113//
114//
115// Initialise the value for P, then determine the actual value depending upon
116// whether the projectile is larger or smaller than the target and these radii
117// in relation to the impact parameter.
118//
119 G4double P = 0.0;
120
121 if (rT > rP)
122 {
123 if (rT-rP<=r && r<=rT+rP) P = 0.125*R*U*S - 0.125*(0.5*R*U+1.0)*T;
124 else P = -1.0;
125 }
126 else
127 {
128 if (rP-rT<=r && r<=rP+rT) P = 0.125*R*U*S - 0.125*(0.5*std::sqrt(n/m)*U-
129 (std::sqrt(1.0-m*m)/n - 1.0)*std::sqrt((2.0-m)/std::pow(m,5.0)))*T;
130 else P = (std::sqrt(1.0-m*m)/n-1.0)*std::sqrt(1.0-b*b/n/n);
131 }
132
133 if (!(P <= 1.0 && P>= -1.0))
134 {
135 if (P > 1.0) P = 1.0;
136 else P = -1.0;
137 }
138 return P;
139}
140////////////////////////////////////////////////////////////////////////////////
141//
142G4double G4NuclearAbrasionGeometry::F ()
143{
144//
145//
146// Initialise the value for F, then determine the actual value depending upon
147// whether the projectile is larger or smaller than the target and these radii
148// in relation to the impact parameter.
149//
150 G4double F = 0.0;
151
152 if (rT > rP)
153 {
154 if (rT-rP<=r && r<=rT+rP) F = 0.75*R*S - 0.125*(3.0*R-1.0)*T;
155 else F = 1.0;
156 }
157 else
158 {
159 if (rP-rT<=r && r<=rP+rT) F = 0.75*R*S - 0.125*(3.0*std::sqrt(n/m)-
160 (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;
161 else F = (1.0-std::pow(1.0-m*m,3.0/2.0))*std::sqrt(1.0-b*b/n/n);
162 }
163
164 if (!(F <= 1.0 && F>= 0.0))
165 {
166 if (F > 1.0) F = 1.0;
167 else F = 0.0;
168 }
169 return F;
170}
171////////////////////////////////////////////////////////////////////////////////
172//
173G4double G4NuclearAbrasionGeometry::GetExcitationEnergyOfProjectile ()
174{
175 G4double F1 = F();
176 G4double P1 = P();
177 G4double Es = 0.0;
178
179 Es = 0.95 * MeV * 4.0 * pi * rP*rP/fermi/fermi *
180 (1.0+P1-std::pow(1.0-F1,2.0/3.0));
181// if (rT < rP && r < rP-rT)
182 if ((r-rP)/rT < rth)
183 {
184 G4double omega = 0.0;
185 if (AP < 12.0) omega = 1500.0;
186 else if (AP <= 16.0) omega = 1500.0 - 320.0*(AP-12.0);
187 Es *= 1.0 + F1*(5.0+omega*F1*F1);
188 }
189
190 if (Es < 0.0)
191 Es = 0.0;
192 else if (Es > B * AP)
193 Es = B * AP;
194 return Es;
195}
196//////////////////////////////////////////////////////////////////////////////////////
197//
198G4double G4NuclearAbrasionGeometry::GetExcitationEnergyOfTarget ()
199{
200//
201//
202// This member function declares a new G4NuclearAbrasionGeometry object but with the
203// projectile and target exchanged to determine the values for F and P. Determination
204// of the excess surface area and excitation energy is as above.
205//
206 G4NuclearAbrasionGeometry *revAbrasionGeometry =
207 new G4NuclearAbrasionGeometry(AT, AP, r);
208 G4double F1 = revAbrasionGeometry->F();
209 G4double P1 = revAbrasionGeometry->P();
210 G4double Es = 0.0;
211
212 Es = 0.95 * MeV * 4.0 * pi * rT*rT/fermi/fermi *
213 (1.0+P1-std::pow(1.0-F1,2.0/3.0));
214// if (rP < rT && r < rT-rP)
215 if ((r-rT)/rP < rth)
216 {
217 G4double omega = 0.0;
218 if (AT < 12.0) omega = 1500.0;
219 else if (AT <= 16.0) omega = 1500.0 - 320.0*(AT-12.0);
220 Es *= 1.0 + F1*(5.0+omega*F1*F1);
221 }
222
223 if (Es < 0.0)
224 Es = 0.0;
225 else if (Es > B * AT)
226 Es = B * AT;
227 return Es;
228}
229////////////////////////////////////////////////////////////////////////////////
230//
Note: See TracBrowser for help on using the repository browser.