ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4NuclearAbrasionGeometry.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4NuclearAbrasionGeometry.cc
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 // 4 August 2015, A. Ribon, CERN, Switzerland
65 // Replacing std::pow with the faster G4Pow.
66 //
67 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
69 //
71 #include "G4WilsonRadius.hh"
72 #include "G4PhysicalConstants.hh"
73 #include "G4SystemOfUnits.hh"
74 #include "G4Pow.hh"
76 //
78  G4double AT1, G4double r1)
79 {
80 //
81 //
82 // Initialise variables for interaction geometry.
83 //
84  G4WilsonRadius aR;
85  AP = AP1;
86  AT = AT1;
87  rP = aR.GetWilsonRadius(AP);
88  rT = aR.GetWilsonRadius(AT);
89  r = r1;
90  n = rP / (rP + rT);
91  b = r / (rP + rT);
92  m = rT / rP;
93  Q = (1.0 - b)/n;
94  S = Q * Q;
95  T = S * Q;
96  R = std::sqrt(m*n);
97  U = 1.0/m - 2.0;
98 //
99 //
100 // Initialise the threshold radius-ratio at which interactions are considered
101 // peripheral or central.
102 //
103  rth = 2.0/3.0;
104  B = 10.0 * MeV;
105 }
107 //
109 {;}
111 //
113  {if (rth1 > 0.0 && rth1 <= 1.0) rth = rth1;}
115 //
117  {return rth;}
119 //
121 {
122 //
123 //
124 // Initialise the value for P, then determine the actual value depending upon
125 // whether the projectile is larger or smaller than the target and these radii
126 // in relation to the impact parameter.
127 //
128  G4double valueP = 0.0;
129 
130  if (rT > rP)
131  {
132  if (rT-rP<=r && r<=rT+rP) valueP = 0.125*R*U*S - 0.125*(0.5*R*U+1.0)*T;
133  else valueP = -1.0;
134  }
135  else
136  {
137  if (rP-rT<=r && r<=rP+rT) valueP = 0.125*R*U*S - 0.125*(0.5*std::sqrt(n/m)*U-
138  (std::sqrt(1.0-m*m)/n - 1.0)*std::sqrt((2.0-m)/G4Pow::GetInstance()->powN(m,5)))*T;
139  else valueP = (std::sqrt(1.0-m*m)/n-1.0)*std::sqrt(1.0-b*b/n/n);
140  }
141 
142  if (!(valueP <= 1.0 && valueP>= -1.0))
143  {
144  if (valueP > 1.0) valueP = 1.0;
145  else valueP = -1.0;
146  }
147  return valueP;
148 }
150 //
152 {
153 //
154 //
155 // Initialise the value for F, then determine the actual value depending upon
156 // whether the projectile is larger or smaller than the target and these radii
157 // in relation to the impact parameter.
158 //
159  G4double valueF = 0.0;
160 
161  if (rT > rP)
162  {
163  if (rT-rP<=r && r<=rT+rP) valueF = 0.75*R*S - 0.125*(3.0*R-1.0)*T;
164  else valueF = 1.0;
165  }
166  else
167  {
168  if (rP-rT<=r && r<=rP+rT) valueF = 0.75*R*S - 0.125*(3.0*std::sqrt(n/m)-
169  (1.0-G4Pow::GetInstance()->powA(1.0-m*m,3.0/2.0))*std::sqrt(1.0-G4Pow::GetInstance()->powN(1.0-m,2))/G4Pow::GetInstance()->powN(m,3))*T;
170  else valueF = (1.0-G4Pow::GetInstance()->powA(1.0-m*m,3.0/2.0))*std::sqrt(1.0-b*b/n/n);
171  }
172 
173  if (!(valueF <= 1.0 && valueF>= 0.0))
174  {
175  if (valueF > 1.0) valueF = 1.0;
176  else valueF = 0.0;
177  }
178  return valueF;
179 }
181 //
183 {
184  G4double F1 = F();
185  G4double P1 = P();
186  G4double Es = 0.0;
187 
188  Es = 0.95 * MeV * 4.0 * pi * rP*rP/fermi/fermi *
189  (1.0+P1-G4Pow::GetInstance()->A23(1.0-F1));
190 // if (rT < rP && r < rP-rT)
191  if ((r-rP)/rT < rth)
192  {
193  G4double omega = 0.0;
194  if (AP < 12.0) omega = 1500.0;
195  else if (AP <= 16.0) omega = 1500.0 - 320.0*(AP-12.0);
196  Es *= 1.0 + F1*(5.0+omega*F1*F1);
197  }
198 
199  if (Es < 0.0)
200  Es = 0.0;
201  else if (Es > B * AP)
202  Es = B * AP;
203  return Es;
204 }
205 
206 
208 {
209  // This member function declares a new G4NuclearAbrasionGeometry object
210  // but with the projectile and target exchanged to determine the values
211  // for F and P. Determination of the excess surface area and excitation
212  // energy is as above.
213 
214  G4NuclearAbrasionGeometry* revAbrasionGeometry =
216  G4double F1 = revAbrasionGeometry->F();
217  G4double P1 = revAbrasionGeometry->P();
218  G4double Es = 0.0;
219 
220  Es = 0.95 * MeV * 4.0 * pi * rT*rT/fermi/fermi *
221  (1.0+P1-G4Pow::GetInstance()->A23(1.0-F1));
222 
223 // if (rP < rT && r < rT-rP)
224  if ((r-rT)/rP < rth) {
225  G4double omega = 0.0;
226  if (AT < 12.0) omega = 1500.0;
227  else if (AT <= 16.0) omega = 1500.0 - 320.0*(AT-12.0);
228  Es *= 1.0 + F1*(5.0+omega*F1*F1);
229  }
230 
231  if (Es < 0.0)
232  Es = 0.0;
233  else if (Es > B * AT)
234  Es = B * AT;
235 
236  delete revAbrasionGeometry;
237 
238  return Es;
239 }