Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. * 20 // * 20 // * * 21 // * Parts of this code which have been devel 21 // * Parts of this code which have been developed by QinetiQ Ltd * 22 // * under contract to the European Space Agen 22 // * under contract to the European Space Agency (ESA) are the * 23 // * intellectual property of ESA. Rights to u 23 // * intellectual property of ESA. Rights to use, copy, modify and * 24 // * redistribute this software for general pu 24 // * redistribute this software for general public use are granted * 25 // * in compliance with any licensing, distrib 25 // * in compliance with any licensing, distribution and development * 26 // * policy adopted by the Geant4 Collaboratio 26 // * policy adopted by the Geant4 Collaboration. This code has been * 27 // * written by QinetiQ Ltd for the European S 27 // * written by QinetiQ Ltd for the European Space Agency, under ESA * 28 // * contract 17191/03/NL/LvH (Aurora Programm 28 // * contract 17191/03/NL/LvH (Aurora Programme). * 29 // * 29 // * * 30 // * By using, copying, modifying or distri 30 // * By using, copying, modifying or distributing the software (or * 31 // * any work based on the software) you ag 31 // * any work based on the software) you agree to acknowledge its * 32 // * use in resulting scientific publicati 32 // * use in resulting scientific publications, and indicate your * 33 // * acceptance of all terms of the Geant4 Sof 33 // * acceptance of all terms of the Geant4 Software license. * 34 // ******************************************* 34 // ******************************************************************** 35 // 35 // 36 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 36 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 37 // 37 // 38 // MODULE: G4NuclearAbrasionGeometry.cc 38 // MODULE: G4NuclearAbrasionGeometry.cc 39 // 39 // 40 // Version: B.1 40 // Version: B.1 41 // Date: 15/04/04 41 // Date: 15/04/04 42 // Author: P R Truscott 42 // Author: P R Truscott 43 // Organisation: QinetiQ Ltd, UK 43 // Organisation: QinetiQ Ltd, UK 44 // Customer: ESA/ESTEC, NOORDWIJK 44 // Customer: ESA/ESTEC, NOORDWIJK 45 // Contract: 17191/03/NL/LvH 45 // Contract: 17191/03/NL/LvH 46 // 46 // 47 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 47 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 48 // 48 // 49 // CHANGE HISTORY 49 // CHANGE HISTORY 50 // -------------- 50 // -------------- 51 // 51 // 52 // 18 November 2003, P R Truscott, QinetiQ Ltd 52 // 18 November 2003, P R Truscott, QinetiQ Ltd, UK 53 // Created. 53 // Created. 54 // 54 // 55 // 15 March 2004, P R Truscott, QinetiQ Ltd, U 55 // 15 March 2004, P R Truscott, QinetiQ Ltd, UK 56 // Beta release 56 // Beta release 57 // 57 // 58 // 4 June 2004, J.P. Wellisch, CERN, Switzerla 58 // 4 June 2004, J.P. Wellisch, CERN, Switzerland 59 // resolving technical portability issues. 59 // resolving technical portability issues. 60 // 60 // 61 // 12 June 2012, A. Ribon, CERN, Switzerland << 62 // Fixing trivial warning errors of shadowed v << 63 // << 64 // 4 August 2015, A. Ribon, CERN, Switzerland << 65 // Replacing std::pow with the faster G4Pow. << 66 // << 67 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 61 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 68 ////////////////////////////////////////////// 62 //////////////////////////////////////////////////////////////////////////////// 69 // 63 // 70 #include "G4NuclearAbrasionGeometry.hh" 64 #include "G4NuclearAbrasionGeometry.hh" 71 #include "G4WilsonRadius.hh" 65 #include "G4WilsonRadius.hh" 72 #include "G4PhysicalConstants.hh" << 73 #include "G4SystemOfUnits.hh" << 74 #include "G4Pow.hh" << 75 ////////////////////////////////////////////// 66 //////////////////////////////////////////////////////////////////////////////// 76 // 67 // 77 G4NuclearAbrasionGeometry::G4NuclearAbrasionGe 68 G4NuclearAbrasionGeometry::G4NuclearAbrasionGeometry (G4double AP1, 78 G4double AT1, G4double r1) 69 G4double AT1, G4double r1) 79 { 70 { 80 // 71 // 81 // 72 // 82 // Initialise variables for interaction geomet 73 // Initialise variables for interaction geometry. 83 // 74 // 84 G4WilsonRadius aR; 75 G4WilsonRadius aR; 85 AP = AP1; 76 AP = AP1; 86 AT = AT1; 77 AT = AT1; 87 rP = aR.GetWilsonRadius(AP); 78 rP = aR.GetWilsonRadius(AP); 88 rT = aR.GetWilsonRadius(AT); 79 rT = aR.GetWilsonRadius(AT); 89 r = r1; 80 r = r1; 90 n = rP / (rP + rT); 81 n = rP / (rP + rT); 91 b = r / (rP + rT); 82 b = r / (rP + rT); 92 m = rT / rP; 83 m = rT / rP; 93 Q = (1.0 - b)/n; 84 Q = (1.0 - b)/n; 94 S = Q * Q; 85 S = Q * Q; 95 T = S * Q; 86 T = S * Q; 96 R = std::sqrt(m*n); 87 R = std::sqrt(m*n); 97 U = 1.0/m - 2.0; 88 U = 1.0/m - 2.0; 98 // 89 // 99 // 90 // 100 // Initialise the threshold radius-ratio at wh 91 // Initialise the threshold radius-ratio at which interactions are considered 101 // peripheral or central. 92 // peripheral or central. 102 // 93 // 103 rth = 2.0/3.0; 94 rth = 2.0/3.0; 104 B = 10.0 * MeV; 95 B = 10.0 * MeV; 105 } 96 } 106 ////////////////////////////////////////////// 97 //////////////////////////////////////////////////////////////////////////////// 107 // 98 // 108 G4NuclearAbrasionGeometry::~G4NuclearAbrasionG 99 G4NuclearAbrasionGeometry::~G4NuclearAbrasionGeometry () 109 {;} 100 {;} 110 ////////////////////////////////////////////// 101 //////////////////////////////////////////////////////////////////////////////// 111 // 102 // 112 void G4NuclearAbrasionGeometry::SetPeripheralT 103 void G4NuclearAbrasionGeometry::SetPeripheralThreshold (G4double rth1) 113 {if (rth1 > 0.0 && rth1 <= 1.0) rth = rth1;} 104 {if (rth1 > 0.0 && rth1 <= 1.0) rth = rth1;} 114 ////////////////////////////////////////////// 105 //////////////////////////////////////////////////////////////////////////////// 115 // 106 // 116 G4double G4NuclearAbrasionGeometry::GetPeriphe 107 G4double G4NuclearAbrasionGeometry::GetPeripheralThreshold () 117 {return rth;} 108 {return rth;} 118 ////////////////////////////////////////////// 109 //////////////////////////////////////////////////////////////////////////////// 119 // 110 // 120 G4double G4NuclearAbrasionGeometry::P () 111 G4double G4NuclearAbrasionGeometry::P () 121 { 112 { 122 // 113 // 123 // 114 // 124 // Initialise the value for P, then determine 115 // Initialise the value for P, then determine the actual value depending upon 125 // whether the projectile is larger or smaller 116 // whether the projectile is larger or smaller than the target and these radii 126 // in relation to the impact parameter. 117 // in relation to the impact parameter. 127 // 118 // 128 G4double valueP = 0.0; << 119 G4double P = 0.0; 129 120 130 if (rT > rP) 121 if (rT > rP) 131 { 122 { 132 if (rT-rP<=r && r<=rT+rP) valueP = 0.125*R << 123 if (rT-rP<=r && r<=rT+rP) P = 0.125*R*U*S - 0.125*(0.5*R*U+1.0)*T; 133 else valueP = -1.0; << 124 else P = -1.0; 134 } 125 } 135 else 126 else 136 { 127 { 137 if (rP-rT<=r && r<=rP+rT) valueP = 0.125*R << 128 if (rP-rT<=r && r<=rP+rT) P = 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(( << 129 (std::sqrt(1.0-m*m)/n - 1.0)*std::sqrt((2.0-m)/std::pow(m,5.0)))*T; 139 else valueP = (std::s << 130 else P = (std::sqrt(1.0-m*m)/n-1.0)*std::sqrt(1.0-b*b/n/n); 140 } 131 } 141 132 142 if (!(valueP <= 1.0 && valueP>= -1.0)) << 133 if (!(P <= 1.0 && P>= -1.0)) 143 { 134 { 144 if (valueP > 1.0) valueP = 1.0; << 135 if (P > 1.0) P = 1.0; 145 else valueP = -1.0; << 136 else P = -1.0; 146 } 137 } 147 return valueP; << 138 return P; 148 } 139 } 149 ////////////////////////////////////////////// 140 //////////////////////////////////////////////////////////////////////////////// 150 // 141 // 151 G4double G4NuclearAbrasionGeometry::F () 142 G4double G4NuclearAbrasionGeometry::F () 152 { 143 { 153 // 144 // 154 // 145 // 155 // Initialise the value for F, then determine 146 // Initialise the value for F, then determine the actual value depending upon 156 // whether the projectile is larger or smaller 147 // whether the projectile is larger or smaller than the target and these radii 157 // in relation to the impact parameter. 148 // in relation to the impact parameter. 158 // 149 // 159 G4double valueF = 0.0; << 150 G4double F = 0.0; 160 151 161 if (rT > rP) 152 if (rT > rP) 162 { 153 { 163 if (rT-rP<=r && r<=rT+rP) valueF = 0.75*R* << 154 if (rT-rP<=r && r<=rT+rP) F = 0.75*R*S - 0.125*(3.0*R-1.0)*T; 164 else valueF = 1.0; << 155 else F = 1.0; 165 } 156 } 166 else 157 else 167 { 158 { 168 if (rP-rT<=r && r<=rP+rT) valueF = 0.75*R* << 159 if (rP-rT<=r && r<=rP+rT) F = 0.75*R*S - 0.125*(3.0*std::sqrt(n/m)- 169 (1.0-G4Pow::GetInstance()->powA(1.0-m*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; 170 else valueF = (1.0-G4 << 161 else F = (1.0-std::pow(1.0-m*m,3.0/2.0))*std::sqrt(1.0-b*b/n/n); 171 } 162 } 172 163 173 if (!(valueF <= 1.0 && valueF>= 0.0)) << 164 if (!(F <= 1.0 && F>= 0.0)) 174 { 165 { 175 if (valueF > 1.0) valueF = 1.0; << 166 if (F > 1.0) F = 1.0; 176 else valueF = 0.0; << 167 else F = 0.0; 177 } 168 } 178 return valueF; << 169 return F; 179 } 170 } 180 ////////////////////////////////////////////// 171 //////////////////////////////////////////////////////////////////////////////// 181 // 172 // 182 G4double G4NuclearAbrasionGeometry::GetExcitat 173 G4double G4NuclearAbrasionGeometry::GetExcitationEnergyOfProjectile () 183 { 174 { 184 G4double F1 = F(); 175 G4double F1 = F(); 185 G4double P1 = P(); 176 G4double P1 = P(); 186 G4double Es = 0.0; 177 G4double Es = 0.0; 187 178 188 Es = 0.95 * MeV * 4.0 * pi * rP*rP/fermi/fer 179 Es = 0.95 * MeV * 4.0 * pi * rP*rP/fermi/fermi * 189 (1.0+P1-G4Pow::GetInstance()->A23(1.0-F << 180 (1.0+P1-std::pow(1.0-F1,2.0/3.0)); 190 // if (rT < rP && r < rP-rT) 181 // if (rT < rP && r < rP-rT) 191 if ((r-rP)/rT < rth) 182 if ((r-rP)/rT < rth) 192 { 183 { 193 G4double omega = 0.0; 184 G4double omega = 0.0; 194 if (AP < 12.0) omega = 1500.0; 185 if (AP < 12.0) omega = 1500.0; 195 else if (AP <= 16.0) omega = 1500.0 - 320. 186 else if (AP <= 16.0) omega = 1500.0 - 320.0*(AP-12.0); 196 Es *= 1.0 + F1*(5.0+omega*F1*F1); 187 Es *= 1.0 + F1*(5.0+omega*F1*F1); 197 } 188 } 198 189 199 if (Es < 0.0) 190 if (Es < 0.0) 200 Es = 0.0; 191 Es = 0.0; 201 else if (Es > B * AP) 192 else if (Es > B * AP) 202 Es = B * AP; 193 Es = B * AP; 203 return Es; 194 return Es; 204 } 195 } 205 << 196 ////////////////////////////////////////////////////////////////////////////////////// 206 << 197 // 207 G4double G4NuclearAbrasionGeometry::GetExcitat 198 G4double G4NuclearAbrasionGeometry::GetExcitationEnergyOfTarget () 208 { 199 { 209 // This member function declares a new G4Nuc << 200 // 210 // but with the projectile and target exchan << 201 // 211 // for F and P. Determination of the excess << 202 // This member function declares a new G4NuclearAbrasionGeometry object but with the 212 // energy is as above. << 203 // projectile and target exchanged to determine the values for F and P. Determination 213 << 204 // of the excess surface area and excitation energy is as above. 214 G4NuclearAbrasionGeometry* revAbrasionGeomet << 205 // >> 206 G4NuclearAbrasionGeometry *revAbrasionGeometry = 215 new G4NuclearAbrasionGeometry(AT, AP, r); 207 new G4NuclearAbrasionGeometry(AT, AP, r); 216 G4double F1 = revAbrasionGeometry->F(); 208 G4double F1 = revAbrasionGeometry->F(); 217 G4double P1 = revAbrasionGeometry->P(); 209 G4double P1 = revAbrasionGeometry->P(); 218 G4double Es = 0.0; 210 G4double Es = 0.0; 219 211 220 Es = 0.95 * MeV * 4.0 * pi * rT*rT/fermi/fer 212 Es = 0.95 * MeV * 4.0 * pi * rT*rT/fermi/fermi * 221 (1.0+P1-G4Pow::GetInstance()->A23(1.0-F << 213 (1.0+P1-std::pow(1.0-F1,2.0/3.0)); 222 << 223 // if (rP < rT && r < rT-rP) 214 // if (rP < rT && r < rT-rP) 224 if ((r-rT)/rP < rth) { << 215 if ((r-rT)/rP < rth) >> 216 { 225 G4double omega = 0.0; 217 G4double omega = 0.0; 226 if (AT < 12.0) omega = 1500.0; 218 if (AT < 12.0) omega = 1500.0; 227 else if (AT <= 16.0) omega = 1500.0 - 320. 219 else if (AT <= 16.0) omega = 1500.0 - 320.0*(AT-12.0); 228 Es *= 1.0 + F1*(5.0+omega*F1*F1); 220 Es *= 1.0 + F1*(5.0+omega*F1*F1); 229 } 221 } 230 222 231 if (Es < 0.0) 223 if (Es < 0.0) 232 Es = 0.0; 224 Es = 0.0; 233 else if (Es > B * AT) 225 else if (Es > B * AT) 234 Es = B * AT; 226 Es = B * AT; 235 << 236 delete revAbrasionGeometry; << 237 << 238 return Es; 227 return Es; 239 } 228 } >> 229 //////////////////////////////////////////////////////////////////////////////// >> 230 // 240 231