| 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: G4WilsonAblationModel. 38 // MODULE: G4WilsonAblationModel.cc
39 // 39 //
40 // Version: 1.0 40 // Version: 1.0
41 // Date: 08/12/2009 41 // Date: 08/12/2009
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 // 6 October 2003, P R Truscott, QinetiQ Ltd, 52 // 6 October 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 // 08 December 2009, P R Truscott, QinetiQ Ltd 58 // 08 December 2009, P R Truscott, QinetiQ Ltd, UK
59 // Ver 1.0 59 // Ver 1.0
60 // Updated as a result of changes in the G4Eva 60 // Updated as a result of changes in the G4Evaporation classes. These changes
61 // affect mostly SelectSecondariesByEvaporatio 61 // affect mostly SelectSecondariesByEvaporation, and now you have variables
62 // associated with the evaporation model which 62 // associated with the evaporation model which can be changed:
63 // OPTxs to select the inverse cross-sectio 63 // OPTxs to select the inverse cross-section
64 // OPTxs = 0 => Dostrovski's parameter 64 // OPTxs = 0 => Dostrovski's parameterization
65 // OPTxs = 1 or 2 => Chatterjee's paramater 65 // OPTxs = 1 or 2 => Chatterjee's paramaterization
66 // OPTxs = 3 or 4 => Kalbach's parameteriza 66 // OPTxs = 3 or 4 => Kalbach's parameterization
67 // useSICB => use superimposed Coulo 67 // useSICB => use superimposed Coulomb Barrier for inverse cross
68 // sections 68 // sections
69 // Other problem found with G4Fragment definit 69 // Other problem found with G4Fragment definition using Lorentz vector and
70 // **G4ParticleDefinition**. This does not al 70 // **G4ParticleDefinition**. This does not allow A and Z to be defined for the
71 // fragment for some reason. Now the fragment 71 // fragment for some reason. Now the fragment is defined more explicitly:
72 // G4Fragment *fragment = new G4Fragment(A, 72 // G4Fragment *fragment = new G4Fragment(A, Z, lorentzVector);
73 // to avoid this quirk. Bug found in SelectSe 73 // to avoid this quirk. Bug found in SelectSecondariesByDefault: *type is now
74 // equated to evapType[i] whereas previously i 74 // equated to evapType[i] whereas previously it was equated to fragType[i].
75 // 75 //
76 // 06 August 2015, A. Ribon, CERN 76 // 06 August 2015, A. Ribon, CERN
77 // Migrated std::exp and std::pow to the faste 77 // Migrated std::exp and std::pow to the faster G4Exp and G4Pow.
78 // 78 //
79 // 09 June 2017, C. Mancini Terracciano, INFN 79 // 09 June 2017, C. Mancini Terracciano, INFN
80 // Fixed bug on the initialization of Photon E 80 // Fixed bug on the initialization of Photon Evaporation model
81 // 81 //
82 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 82 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
83 ////////////////////////////////////////////// 83 ////////////////////////////////////////////////////////////////////////////////
84 // 84 //
85 #include <iomanip> 85 #include <iomanip>
86 #include <numeric> 86 #include <numeric>
87 87
88 #include "G4WilsonAblationModel.hh" 88 #include "G4WilsonAblationModel.hh"
89 #include "G4PhysicalConstants.hh" 89 #include "G4PhysicalConstants.hh"
90 #include "G4SystemOfUnits.hh" 90 #include "G4SystemOfUnits.hh"
91 #include "Randomize.hh" 91 #include "Randomize.hh"
92 #include "G4ParticleTable.hh" 92 #include "G4ParticleTable.hh"
93 #include "G4IonTable.hh" 93 #include "G4IonTable.hh"
94 #include "G4Alpha.hh" 94 #include "G4Alpha.hh"
95 #include "G4He3.hh" 95 #include "G4He3.hh"
96 #include "G4Triton.hh" 96 #include "G4Triton.hh"
97 #include "G4Deuteron.hh" 97 #include "G4Deuteron.hh"
98 #include "G4Proton.hh" 98 #include "G4Proton.hh"
99 #include "G4Neutron.hh" 99 #include "G4Neutron.hh"
100 #include "G4AlphaEvaporationChannel.hh" 100 #include "G4AlphaEvaporationChannel.hh"
101 #include "G4He3EvaporationChannel.hh" 101 #include "G4He3EvaporationChannel.hh"
102 #include "G4TritonEvaporationChannel.hh" 102 #include "G4TritonEvaporationChannel.hh"
103 #include "G4DeuteronEvaporationChannel.hh" 103 #include "G4DeuteronEvaporationChannel.hh"
104 #include "G4ProtonEvaporationChannel.hh" 104 #include "G4ProtonEvaporationChannel.hh"
105 #include "G4NeutronEvaporationChannel.hh" 105 #include "G4NeutronEvaporationChannel.hh"
106 #include "G4PhotonEvaporation.hh" 106 #include "G4PhotonEvaporation.hh"
107 #include "G4LorentzVector.hh" 107 #include "G4LorentzVector.hh"
108 #include "G4VEvaporationChannel.hh" 108 #include "G4VEvaporationChannel.hh"
109 109
110 #include "G4Exp.hh" 110 #include "G4Exp.hh"
111 #include "G4Pow.hh" 111 #include "G4Pow.hh"
112 112
113 #include "G4PhysicsModelCatalog.hh" <<
114 113
115 ////////////////////////////////////////////// 114 ////////////////////////////////////////////////////////////////////////////////
116 // 115 //
117 G4WilsonAblationModel::G4WilsonAblationModel() 116 G4WilsonAblationModel::G4WilsonAblationModel()
118 { 117 {
119 // 118 //
120 // 119 //
121 // Send message to stdout to advise that the G 120 // Send message to stdout to advise that the G4Abrasion model is being used.
122 // 121 //
123 PrintWelcomeMessage(); 122 PrintWelcomeMessage();
124 // 123 //
125 // 124 //
126 // Set the default verbose level to 0 - no out 125 // Set the default verbose level to 0 - no output.
127 // 126 //
128 verboseLevel = 0; 127 verboseLevel = 0;
129 // 128 //
130 // 129 //
131 // Set the binding energy per nucleon .... did 130 // Set the binding energy per nucleon .... did I mention that this is a crude
132 // model for nuclear de-excitation? 131 // model for nuclear de-excitation?
133 // 132 //
134 B = 10.0 * MeV; 133 B = 10.0 * MeV;
135 // 134 //
136 // 135 //
137 // It is possuble to switch off secondary part 136 // It is possuble to switch off secondary particle production (other than the
138 // final nuclear fragment). The default is on 137 // final nuclear fragment). The default is on.
139 // 138 //
140 produceSecondaries = true; 139 produceSecondaries = true;
141 // 140 //
142 // 141 //
143 // Now we need to define the decay modes. We' 142 // Now we need to define the decay modes. We're using the G4Evaporation model
144 // to help determine the kinematics of the dec 143 // to help determine the kinematics of the decay.
145 // 144 //
146 nFragTypes = 6; 145 nFragTypes = 6;
147 fragType[0] = G4Alpha::Alpha(); 146 fragType[0] = G4Alpha::Alpha();
148 fragType[1] = G4He3::He3(); 147 fragType[1] = G4He3::He3();
149 fragType[2] = G4Triton::Triton(); 148 fragType[2] = G4Triton::Triton();
150 fragType[3] = G4Deuteron::Deuteron(); 149 fragType[3] = G4Deuteron::Deuteron();
151 fragType[4] = G4Proton::Proton(); 150 fragType[4] = G4Proton::Proton();
152 fragType[5] = G4Neutron::Neutron(); 151 fragType[5] = G4Neutron::Neutron();
153 for(G4int i=0; i<200; ++i) { fSig[i] = 0.0; 152 for(G4int i=0; i<200; ++i) { fSig[i] = 0.0; }
154 // 153 //
155 // 154 //
156 // Set verboseLevel default to no output. 155 // Set verboseLevel default to no output.
157 // 156 //
158 verboseLevel = 0; 157 verboseLevel = 0;
159 theChannelFactory = new G4EvaporationFactory 158 theChannelFactory = new G4EvaporationFactory(new G4PhotonEvaporation());
160 theChannels = theChannelFactory->GetChannel( 159 theChannels = theChannelFactory->GetChannel();
161 // 160 //
162 // 161 //
163 // Set defaults for evaporation classes. Thes 162 // Set defaults for evaporation classes. These can be overridden by user
164 // "set" methods. 163 // "set" methods.
165 // 164 //
166 OPTxs = 3; 165 OPTxs = 3;
167 useSICB = false; 166 useSICB = false;
168 fragmentVector = 0; 167 fragmentVector = 0;
169 <<
170 secID = G4PhysicsModelCatalog::GetModelID("m <<
171 } 168 }
172 ////////////////////////////////////////////// 169 ////////////////////////////////////////////////////////////////////////////////
173 // 170 //
174 G4WilsonAblationModel::~G4WilsonAblationModel( 171 G4WilsonAblationModel::~G4WilsonAblationModel()
175 {} 172 {}
176 173
177 ////////////////////////////////////////////// 174 ////////////////////////////////////////////////////////////////////////////////
178 // 175 //
179 G4FragmentVector *G4WilsonAblationModel::Break 176 G4FragmentVector *G4WilsonAblationModel::BreakItUp
180 (const G4Fragment &theNucleus) 177 (const G4Fragment &theNucleus)
181 { 178 {
182 // 179 //
183 // 180 //
184 // Initilise the pointer to the G4FragmentVect 181 // Initilise the pointer to the G4FragmentVector used to return the information
185 // about the breakup. 182 // about the breakup.
186 // 183 //
187 fragmentVector = new G4FragmentVector; 184 fragmentVector = new G4FragmentVector;
188 fragmentVector->clear(); 185 fragmentVector->clear();
189 // 186 //
190 // 187 //
191 // Get the A, Z and excitation of the nucleus. 188 // Get the A, Z and excitation of the nucleus.
192 // 189 //
193 G4int A = theNucleus.GetA_asInt(); 190 G4int A = theNucleus.GetA_asInt();
194 G4int Z = theNucleus.GetZ_asInt(); 191 G4int Z = theNucleus.GetZ_asInt();
195 G4double ex = theNucleus.GetExcitationEnergy 192 G4double ex = theNucleus.GetExcitationEnergy();
196 if (verboseLevel >= 2) 193 if (verboseLevel >= 2)
197 { 194 {
198 G4cout <<"oooooooooooooooooooooooooooooooo 195 G4cout <<"oooooooooooooooooooooooooooooooooooooooo"
199 <<"oooooooooooooooooooooooooooooooo 196 <<"oooooooooooooooooooooooooooooooooooooooo"
200 <<G4endl; 197 <<G4endl;
201 G4cout.precision(6); 198 G4cout.precision(6);
202 G4cout <<"IN G4WilsonAblationModel" <<G4en 199 G4cout <<"IN G4WilsonAblationModel" <<G4endl;
203 G4cout <<"Initial prefragment A=" <<A 200 G4cout <<"Initial prefragment A=" <<A
204 <<", Z=" <<Z 201 <<", Z=" <<Z
205 <<", excitation energy = " <<ex/MeV 202 <<", excitation energy = " <<ex/MeV <<" MeV"
206 <<G4endl; 203 <<G4endl;
207 } 204 }
208 // 205 //
209 // 206 //
210 // Check that there is a nucleus to speak of. 207 // Check that there is a nucleus to speak of. It's possible there isn't one
211 // or its just a proton or neutron. In either 208 // or its just a proton or neutron. In either case, the excitation energy
212 // (from the Lorentz vector) is not used. 209 // (from the Lorentz vector) is not used.
213 // 210 //
214 if (A == 0) 211 if (A == 0)
215 { 212 {
216 if (verboseLevel >= 2) 213 if (verboseLevel >= 2)
217 { 214 {
218 G4cout <<"No nucleus to decay" <<G4endl; 215 G4cout <<"No nucleus to decay" <<G4endl;
219 G4cout <<"oooooooooooooooooooooooooooooo 216 G4cout <<"oooooooooooooooooooooooooooooooooooooooo"
220 <<"oooooooooooooooooooooooooooooo 217 <<"oooooooooooooooooooooooooooooooooooooooo"
221 <<G4endl; 218 <<G4endl;
222 } 219 }
223 return fragmentVector; 220 return fragmentVector;
224 } 221 }
225 else if (A == 1) 222 else if (A == 1)
226 { 223 {
227 G4LorentzVector lorentzVector = theNucleus 224 G4LorentzVector lorentzVector = theNucleus.GetMomentum();
228 lorentzVector.setE(lorentzVector.e()-ex+10 225 lorentzVector.setE(lorentzVector.e()-ex+10.0*eV);
229 if (Z == 0) 226 if (Z == 0)
230 { 227 {
231 G4Fragment *fragment = new G4Fragment(lo 228 G4Fragment *fragment = new G4Fragment(lorentzVector,G4Neutron::Neutron());
232 if (fragment != nullptr) { fragment->Set <<
233 fragmentVector->push_back(fragment); 229 fragmentVector->push_back(fragment);
234 } 230 }
235 else 231 else
236 { 232 {
237 G4Fragment *fragment = new G4Fragment(lo 233 G4Fragment *fragment = new G4Fragment(lorentzVector,G4Proton::Proton());
238 if (fragment != nullptr) { fragment->Set <<
239 fragmentVector->push_back(fragment); 234 fragmentVector->push_back(fragment);
240 } 235 }
241 if (verboseLevel >= 2) 236 if (verboseLevel >= 2)
242 { 237 {
243 G4cout <<"Final fragment is in fact only 238 G4cout <<"Final fragment is in fact only a nucleon) :" <<G4endl;
244 G4cout <<(*fragmentVector)[0] <<G4endl; 239 G4cout <<(*fragmentVector)[0] <<G4endl;
245 G4cout <<"oooooooooooooooooooooooooooooo 240 G4cout <<"oooooooooooooooooooooooooooooooooooooooo"
246 <<"oooooooooooooooooooooooooooooo 241 <<"oooooooooooooooooooooooooooooooooooooooo"
247 <<G4endl; 242 <<G4endl;
248 } 243 }
249 return fragmentVector; 244 return fragmentVector;
250 } 245 }
251 // 246 //
252 // 247 //
253 // Then the number of nucleons ablated (either 248 // Then the number of nucleons ablated (either as nucleons or light nuclear
254 // fragments) is based on a simple argument fo 249 // fragments) is based on a simple argument for the binding energy per nucleon.
255 // 250 //
256 G4int DAabl = (G4int) (ex / B); 251 G4int DAabl = (G4int) (ex / B);
257 if (DAabl > A) DAabl = A; 252 if (DAabl > A) DAabl = A;
258 // The following lines are no longer accurate 253 // The following lines are no longer accurate given we now treat the final fragment
259 // if (verboseLevel >= 2) 254 // if (verboseLevel >= 2)
260 // G4cout <<"Number of nucleons ejected = " 255 // G4cout <<"Number of nucleons ejected = " <<DAabl <<G4endl;
261 256
262 // 257 //
263 // 258 //
264 // Determine the nuclear fragment from the abl 259 // Determine the nuclear fragment from the ablation process by sampling the
265 // Rudstam equation. 260 // Rudstam equation.
266 // 261 //
267 G4int AF = A - DAabl; 262 G4int AF = A - DAabl;
268 G4int ZF = 0; 263 G4int ZF = 0;
269 264
270 if (AF > 0) 265 if (AF > 0)
271 { 266 {
272 G4Pow* g4calc = G4Pow::GetInstance(); 267 G4Pow* g4calc = G4Pow::GetInstance();
273 G4double AFd = (G4double) AF; 268 G4double AFd = (G4double) AF;
274 G4double R = 11.8 / g4calc->powZ(AF, 0.45) 269 G4double R = 11.8 / g4calc->powZ(AF, 0.45);
275 G4int minZ = std::max(1, Z - DAabl); 270 G4int minZ = std::max(1, Z - DAabl);
276 // 271 //
277 // 272 //
278 // Here we define an integral probability dist 273 // Here we define an integral probability distribution based on the Rudstam
279 // equation assuming a constant AF. 274 // equation assuming a constant AF.
280 // 275 //
281 G4int zmax = std::min(199, Z); 276 G4int zmax = std::min(199, Z);
282 G4double sum = 0.0; 277 G4double sum = 0.0;
283 for (ZF=minZ; ZF<=zmax; ++ZF) 278 for (ZF=minZ; ZF<=zmax; ++ZF)
284 { 279 {
285 sum += G4Exp(-R*g4calc->powA(std::abs(ZF 280 sum += G4Exp(-R*g4calc->powA(std::abs(ZF - 0.486*AFd + 3.8E-04*AFd*AFd),1.5));
286 fSig[ZF] = sum; 281 fSig[ZF] = sum;
287 } 282 }
288 // 283 //
289 // 284 //
290 // Now sample that distribution to determine a 285 // Now sample that distribution to determine a value for ZF.
291 // 286 //
292 sum *= G4UniformRand(); 287 sum *= G4UniformRand();
293 for (ZF=minZ; ZF<=zmax; ++ZF) { 288 for (ZF=minZ; ZF<=zmax; ++ZF) {
294 if(sum <= fSig[ZF]) { break; } 289 if(sum <= fSig[ZF]) { break; }
295 } 290 }
296 } 291 }
297 G4int DZabl = Z - ZF; 292 G4int DZabl = Z - ZF;
298 // 293 //
299 // 294 //
300 // Now determine the nucleons or nuclei which 295 // Now determine the nucleons or nuclei which have bee ablated. The preference
301 // is for the production of alphas, then other 296 // is for the production of alphas, then other nuclei in order of decreasing
302 // binding energy. The energies assigned to th 297 // binding energy. The energies assigned to the products of the decay are
303 // provisional for the moment (the 10eV is jus 298 // provisional for the moment (the 10eV is just to avoid errors with negative
304 // excitation energies due to rounding). 299 // excitation energies due to rounding).
305 // 300 //
306 G4double totalEpost = 0.0; 301 G4double totalEpost = 0.0;
307 evapType.clear(); 302 evapType.clear();
308 for (G4int ift=0; ift<nFragTypes; ift++) 303 for (G4int ift=0; ift<nFragTypes; ift++)
309 { 304 {
310 G4ParticleDefinition *type = fragType[ift] 305 G4ParticleDefinition *type = fragType[ift];
311 G4double n = std::floor((G4double) DAabl 306 G4double n = std::floor((G4double) DAabl / type->GetBaryonNumber() + 1.0E-10);
312 G4double n1 = 1.0E+10; 307 G4double n1 = 1.0E+10;
313 if (fragType[ift]->GetPDGCharge() > 0.0) 308 if (fragType[ift]->GetPDGCharge() > 0.0)
314 n1 = std::floor((G4double) DZabl / type- 309 n1 = std::floor((G4double) DZabl / type->GetPDGCharge() + 1.0E-10);
315 if (n > n1) n = n1; 310 if (n > n1) n = n1;
316 if (n > 0.0) 311 if (n > 0.0)
317 { 312 {
318 G4double mass = type->GetPDGMass(); 313 G4double mass = type->GetPDGMass();
319 for (G4int j=0; j<(G4int) n; j++) 314 for (G4int j=0; j<(G4int) n; j++)
320 { 315 {
321 totalEpost += mass; 316 totalEpost += mass;
322 evapType.push_back(type); 317 evapType.push_back(type);
323 } 318 }
324 DAabl -= (G4int) (n * type->GetBaryonNum 319 DAabl -= (G4int) (n * type->GetBaryonNumber() + 1.0E-10);
325 DZabl -= (G4int) (n * type->GetPDGCharge 320 DZabl -= (G4int) (n * type->GetPDGCharge() + 1.0E-10);
326 } 321 }
327 } 322 }
328 // 323 //
329 // 324 //
330 // Determine the properties of the final nucle 325 // Determine the properties of the final nuclear fragment. Note that if
331 // the final fragment is predicted to have a n 326 // the final fragment is predicted to have a nucleon number of zero, then
332 // really it's the particle last in the vector 327 // really it's the particle last in the vector evapType which becomes the
333 // final fragment. Therefore delete this from 328 // final fragment. Therefore delete this from the vector if this is the
334 // case. 329 // case.
335 // 330 //
336 G4double massFinalFrag = 0.0; 331 G4double massFinalFrag = 0.0;
337 if (AF > 0) 332 if (AF > 0)
338 massFinalFrag = G4ParticleTable::GetPartic 333 massFinalFrag = G4ParticleTable::GetParticleTable()->GetIonTable()->
339 GetIonMass(ZF,AF); 334 GetIonMass(ZF,AF);
340 else 335 else
341 { 336 {
342 G4ParticleDefinition *type = evapType[evap 337 G4ParticleDefinition *type = evapType[evapType.size()-1];
343 AF = type->GetBary 338 AF = type->GetBaryonNumber();
344 ZF = (G4int) (type 339 ZF = (G4int) (type->GetPDGCharge() + 1.0E-10);
345 evapType.erase(evapType.end()-1); 340 evapType.erase(evapType.end()-1);
346 } 341 }
347 totalEpost += massFinalFrag; 342 totalEpost += massFinalFrag;
348 // 343 //
349 // 344 //
350 // Provide verbose output on the nuclear fragm 345 // Provide verbose output on the nuclear fragment if requested.
351 // 346 //
352 if (verboseLevel >= 2) 347 if (verboseLevel >= 2)
353 { 348 {
354 G4cout <<"Final fragment A=" <<AF 349 G4cout <<"Final fragment A=" <<AF
355 <<", Z=" <<ZF 350 <<", Z=" <<ZF
356 <<G4endl; 351 <<G4endl;
357 for (G4int ift=0; ift<nFragTypes; ift++) 352 for (G4int ift=0; ift<nFragTypes; ift++)
358 { 353 {
359 G4ParticleDefinition *type = fragType[if 354 G4ParticleDefinition *type = fragType[ift];
360 G4long n = std::count(evapType.cbegin(), << 355 G4int n = std::count(evapType.begin(),evapType.end(),type);
361 if (n > 0) 356 if (n > 0)
362 G4cout <<"Particle type: " <<std::setw 357 G4cout <<"Particle type: " <<std::setw(10) <<type->GetParticleName()
363 <<", number of particles emitte 358 <<", number of particles emitted = " <<n <<G4endl;
364 } 359 }
365 } 360 }
366 // 361 //
367 // Add the total energy from the fragment. No 362 // Add the total energy from the fragment. Note that the fragment is assumed
368 // to be de-excited and does not undergo photo 363 // to be de-excited and does not undergo photo-evaporation .... I did mention
369 // this is a bit of a crude model? 364 // this is a bit of a crude model?
370 // 365 //
371 G4double massPreFrag = theNucleus.GetGr 366 G4double massPreFrag = theNucleus.GetGroundStateMass();
372 G4double totalEpre = massPreFrag + ex 367 G4double totalEpre = massPreFrag + ex;
373 G4double excess = totalEpre - tota 368 G4double excess = totalEpre - totalEpost;
374 // G4Fragment *resultNucleus(theNucleus); 369 // G4Fragment *resultNucleus(theNucleus);
375 G4Fragment *resultNucleus = new G4Fragment(A 370 G4Fragment *resultNucleus = new G4Fragment(A, Z, theNucleus.GetMomentum());
376 G4ThreeVector boost(0.0,0.0,0.0); 371 G4ThreeVector boost(0.0,0.0,0.0);
377 std::size_t nEvap = 0; << 372 G4int nEvap = 0;
378 if (produceSecondaries && evapType.size()>0) 373 if (produceSecondaries && evapType.size()>0)
379 { 374 {
380 if (excess > 0.0) 375 if (excess > 0.0)
381 { 376 {
382 SelectSecondariesByEvaporation (resultNu 377 SelectSecondariesByEvaporation (resultNucleus);
383 nEvap = fragmentVector->size(); 378 nEvap = fragmentVector->size();
384 boost = resultNucleus->GetMomentum().fin 379 boost = resultNucleus->GetMomentum().findBoostToCM();
385 if (evapType.size() > 0) 380 if (evapType.size() > 0)
386 SelectSecondariesByDefault (boost); 381 SelectSecondariesByDefault (boost);
387 } 382 }
388 else 383 else
389 SelectSecondariesByDefault(G4ThreeVector 384 SelectSecondariesByDefault(G4ThreeVector(0.0,0.0,0.0));
390 } 385 }
391 386
392 if (AF > 0) 387 if (AF > 0)
393 { 388 {
394 G4double mass = G4ParticleTable::GetPartic 389 G4double mass = G4ParticleTable::GetParticleTable()->GetIonTable()->
395 GetIonMass(ZF,AF); 390 GetIonMass(ZF,AF);
396 G4double e = mass + 10.0*eV; 391 G4double e = mass + 10.0*eV;
397 G4double p = std::sqrt(e*e-mass*mass); 392 G4double p = std::sqrt(e*e-mass*mass);
398 G4ThreeVector direction(0.0,0.0,1.0); 393 G4ThreeVector direction(0.0,0.0,1.0);
399 G4LorentzVector lorentzVector = G4LorentzV 394 G4LorentzVector lorentzVector = G4LorentzVector(direction*p, e);
400 lorentzVector.boost(-boost); 395 lorentzVector.boost(-boost);
401 G4Fragment* frag = new G4Fragment(AF, ZF, 396 G4Fragment* frag = new G4Fragment(AF, ZF, lorentzVector);
402 if (frag != nullptr) { frag->SetCreatorMod <<
403 fragmentVector->push_back(frag); 397 fragmentVector->push_back(frag);
404 } 398 }
405 delete resultNucleus; 399 delete resultNucleus;
406 // 400 //
407 // 401 //
408 // Provide verbose output on the ablation prod 402 // Provide verbose output on the ablation products if requested.
409 // 403 //
410 if (verboseLevel >= 2) 404 if (verboseLevel >= 2)
411 { 405 {
412 if (nEvap > 0) 406 if (nEvap > 0)
413 { 407 {
414 G4cout <<"----------------------" <<G4en 408 G4cout <<"----------------------" <<G4endl;
415 G4cout <<"Evaporated particles :" <<G4en 409 G4cout <<"Evaporated particles :" <<G4endl;
416 G4cout <<"----------------------" <<G4en 410 G4cout <<"----------------------" <<G4endl;
417 } 411 }
418 std::size_t ie = 0; << 412 G4int ie = 0;
419 for (auto iter = fragmentVector->cbegin(); << 413 G4FragmentVector::iterator iter;
420 iter != fragmentVector->cend(); << 414 for (iter = fragmentVector->begin(); iter != fragmentVector->end(); iter++)
421 { 415 {
422 if (ie == nEvap) 416 if (ie == nEvap)
423 { 417 {
424 // G4cout <<*iter <<G4endl; 418 // G4cout <<*iter <<G4endl;
425 G4cout <<"---------------------------- 419 G4cout <<"---------------------------------" <<G4endl;
426 G4cout <<"Particles from default emiss 420 G4cout <<"Particles from default emission :" <<G4endl;
427 G4cout <<"---------------------------- 421 G4cout <<"---------------------------------" <<G4endl;
428 } 422 }
429 G4cout <<*iter <<G4endl; 423 G4cout <<*iter <<G4endl;
430 } 424 }
431 G4cout <<"oooooooooooooooooooooooooooooooo 425 G4cout <<"oooooooooooooooooooooooooooooooooooooooo"
432 <<"oooooooooooooooooooooooooooooooo 426 <<"oooooooooooooooooooooooooooooooooooooooo"
433 <<G4endl; 427 <<G4endl;
434 } 428 }
435 429
436 return fragmentVector; 430 return fragmentVector;
437 } 431 }
438 ////////////////////////////////////////////// 432 ////////////////////////////////////////////////////////////////////////////////
439 // 433 //
440 void G4WilsonAblationModel::SelectSecondariesB 434 void G4WilsonAblationModel::SelectSecondariesByEvaporation
441 (G4Fragment *intermediateNucleus) 435 (G4Fragment *intermediateNucleus)
442 { 436 {
443 G4Fragment theResidualNucleus = *intermediat 437 G4Fragment theResidualNucleus = *intermediateNucleus;
444 G4bool evaporate = true; 438 G4bool evaporate = true;
445 // Loop checking, 05-Aug-2015, Vladimir Ivan 439 // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
446 while (evaporate && evapType.size() != 0) 440 while (evaporate && evapType.size() != 0)
447 { 441 {
448 // 442 //
449 // 443 //
450 // Here's the cheaky bit. We're hijacking the 444 // Here's the cheaky bit. We're hijacking the G4Evaporation model, in order to
451 // more accurately sample to kinematics, but t 445 // more accurately sample to kinematics, but the species of the nuclear
452 // fragments will be the ones of our choosing 446 // fragments will be the ones of our choosing as above.
453 // 447 //
454 std::vector <G4VEvaporationChannel*> theC 448 std::vector <G4VEvaporationChannel*> theChannels1;
455 theChannels1.clear(); 449 theChannels1.clear();
456 std::vector <G4VEvaporationChannel*>::iter 450 std::vector <G4VEvaporationChannel*>::iterator i;
457 VectorOfFragmentTypes::iterator iter; 451 VectorOfFragmentTypes::iterator iter;
458 std::vector <VectorOfFragmentTypes::iterat 452 std::vector <VectorOfFragmentTypes::iterator> iters;
459 iters.clear(); 453 iters.clear();
460 iter = std::find(evapType.begin(), evapTyp 454 iter = std::find(evapType.begin(), evapType.end(), G4Alpha::Alpha());
461 if (iter != evapType.end()) 455 if (iter != evapType.end())
462 { 456 {
463 theChannels1.push_back(new G4AlphaEvapor 457 theChannels1.push_back(new G4AlphaEvaporationChannel);
464 i = theChannels1.end() - 1; 458 i = theChannels1.end() - 1;
465 (*i)->SetOPTxs(OPTxs); 459 (*i)->SetOPTxs(OPTxs);
466 (*i)->UseSICB(useSICB); 460 (*i)->UseSICB(useSICB);
467 // (*i)->Initialize(theResidualNucleus); 461 // (*i)->Initialize(theResidualNucleus);
468 iters.push_back(iter); 462 iters.push_back(iter);
469 } 463 }
470 iter = std::find(evapType.begin(), evapTyp 464 iter = std::find(evapType.begin(), evapType.end(), G4He3::He3());
471 if (iter != evapType.end()) 465 if (iter != evapType.end())
472 { 466 {
473 theChannels1.push_back(new G4He3Evaporat 467 theChannels1.push_back(new G4He3EvaporationChannel);
474 i = theChannels1.end() - 1; 468 i = theChannels1.end() - 1;
475 (*i)->SetOPTxs(OPTxs); 469 (*i)->SetOPTxs(OPTxs);
476 (*i)->UseSICB(useSICB); 470 (*i)->UseSICB(useSICB);
477 // (*i)->Initialize(theResidualNucleus); 471 // (*i)->Initialize(theResidualNucleus);
478 iters.push_back(iter); 472 iters.push_back(iter);
479 } 473 }
480 iter = std::find(evapType.begin(), evapTyp 474 iter = std::find(evapType.begin(), evapType.end(), G4Triton::Triton());
481 if (iter != evapType.end()) 475 if (iter != evapType.end())
482 { 476 {
483 theChannels1.push_back(new G4TritonEvapo 477 theChannels1.push_back(new G4TritonEvaporationChannel);
484 i = theChannels1.end() - 1; 478 i = theChannels1.end() - 1;
485 (*i)->SetOPTxs(OPTxs); 479 (*i)->SetOPTxs(OPTxs);
486 (*i)->UseSICB(useSICB); 480 (*i)->UseSICB(useSICB);
487 // (*i)->Initialize(theResidualNucleus); 481 // (*i)->Initialize(theResidualNucleus);
488 iters.push_back(iter); 482 iters.push_back(iter);
489 } 483 }
490 iter = std::find(evapType.begin(), evapTyp 484 iter = std::find(evapType.begin(), evapType.end(), G4Deuteron::Deuteron());
491 if (iter != evapType.end()) 485 if (iter != evapType.end())
492 { 486 {
493 theChannels1.push_back(new G4DeuteronEva 487 theChannels1.push_back(new G4DeuteronEvaporationChannel);
494 i = theChannels1.end() - 1; 488 i = theChannels1.end() - 1;
495 (*i)->SetOPTxs(OPTxs); 489 (*i)->SetOPTxs(OPTxs);
496 (*i)->UseSICB(useSICB); 490 (*i)->UseSICB(useSICB);
497 // (*i)->Initialize(theResidualNucleus); 491 // (*i)->Initialize(theResidualNucleus);
498 iters.push_back(iter); 492 iters.push_back(iter);
499 } 493 }
500 iter = std::find(evapType.begin(), evapTyp 494 iter = std::find(evapType.begin(), evapType.end(), G4Proton::Proton());
501 if (iter != evapType.end()) 495 if (iter != evapType.end())
502 { 496 {
503 theChannels1.push_back(new G4ProtonEvapo 497 theChannels1.push_back(new G4ProtonEvaporationChannel);
504 i = theChannels1.end() - 1; 498 i = theChannels1.end() - 1;
505 (*i)->SetOPTxs(OPTxs); 499 (*i)->SetOPTxs(OPTxs);
506 (*i)->UseSICB(useSICB); 500 (*i)->UseSICB(useSICB);
507 // (*i)->Initialize(theResidualNucleus); 501 // (*i)->Initialize(theResidualNucleus);
508 iters.push_back(iter); 502 iters.push_back(iter);
509 } 503 }
510 iter = std::find(evapType.begin(), evapTyp 504 iter = std::find(evapType.begin(), evapType.end(), G4Neutron::Neutron());
511 if (iter != evapType.end()) 505 if (iter != evapType.end())
512 { 506 {
513 theChannels1.push_back(new G4NeutronEvap 507 theChannels1.push_back(new G4NeutronEvaporationChannel);
514 i = theChannels1.end() - 1; 508 i = theChannels1.end() - 1;
515 (*i)->SetOPTxs(OPTxs); 509 (*i)->SetOPTxs(OPTxs);
516 (*i)->UseSICB(useSICB); 510 (*i)->UseSICB(useSICB);
517 // (*i)->Initialize(theResidualNucleus); 511 // (*i)->Initialize(theResidualNucleus);
518 iters.push_back(iter); 512 iters.push_back(iter);
519 } 513 }
520 std::size_t nChannels = theChannels1.size( << 514 G4int nChannels = theChannels1.size();
521 515
522 G4double totalProb = 0.0; 516 G4double totalProb = 0.0;
523 G4int ich = 0; 517 G4int ich = 0;
524 G4double probEvapType[6] = {0.0}; 518 G4double probEvapType[6] = {0.0};
525 for (auto iterEv=theChannels1.cbegin(); << 519 std::vector<G4VEvaporationChannel*>::iterator iterEv;
526 iterEv!=theChannels1.cend(); ++i << 520 for (iterEv=theChannels1.begin(); iterEv!=theChannels1.end(); iterEv++) {
527 totalProb += (*iterEv)->GetEmissionProba 521 totalProb += (*iterEv)->GetEmissionProbability(intermediateNucleus);
528 probEvapType[ich] = totalProb; 522 probEvapType[ich] = totalProb;
529 ++ich; 523 ++ich;
530 } 524 }
531 if (totalProb > 0.0) { 525 if (totalProb > 0.0) {
532 // 526 //
533 // 527 //
534 // The emission probability for at least one o 528 // The emission probability for at least one of the evaporation channels is
535 // positive, therefore work out which one shou 529 // positive, therefore work out which one should be selected and decay
536 // the nucleus. 530 // the nucleus.
537 // 531 //
538 G4double xi = totalProb*G4UniformRand(); 532 G4double xi = totalProb*G4UniformRand();
539 std::size_t ii = 0; << 533 G4int ii = 0;
540 for (ii=0; ii<nChannels; ++ii) << 534 for (ii=0; ii<nChannels; ii++) {
541 { <<
542 if (xi < probEvapType[ii]) { break; } 535 if (xi < probEvapType[ii]) { break; }
543 } 536 }
544 if (ii >= nChannels) { ii = nChannels - 537 if (ii >= nChannels) { ii = nChannels - 1; }
545 G4FragmentVector *evaporationResult = th 538 G4FragmentVector *evaporationResult = theChannels1[ii]->
546 BreakUpFragment(intermediateNucleus); 539 BreakUpFragment(intermediateNucleus);
547 if ((*evaporationResult)[0] != nullptr) <<
548 { <<
549 (*evaporationResult)[0]->SetCreatorMod <<
550 } <<
551 fragmentVector->push_back((*evaporationR 540 fragmentVector->push_back((*evaporationResult)[0]);
552 intermediateNucleus = (*evaporationResul 541 intermediateNucleus = (*evaporationResult)[1];
553 delete evaporationResult; 542 delete evaporationResult;
554 } 543 }
555 else 544 else
556 { 545 {
557 // 546 //
558 // 547 //
559 // Probability for further evaporation is nil 548 // Probability for further evaporation is nil so have to escape from this
560 // routine and set the energies of the seconda 549 // routine and set the energies of the secondaries to 10eV.
561 // 550 //
562 evaporate = false; 551 evaporate = false;
563 } 552 }
564 } 553 }
565 554
566 return; 555 return;
567 } 556 }
568 ////////////////////////////////////////////// 557 ////////////////////////////////////////////////////////////////////////////////
569 // 558 //
570 void G4WilsonAblationModel::SelectSecondariesB 559 void G4WilsonAblationModel::SelectSecondariesByDefault (G4ThreeVector boost)
571 { 560 {
572 for (std::size_t i=0; i<evapType.size(); ++i << 561 for (unsigned i=0; i<evapType.size(); i++)
573 { 562 {
574 G4ParticleDefinition *type = evapType[i]; 563 G4ParticleDefinition *type = evapType[i];
575 G4double mass = type->GetPDGM 564 G4double mass = type->GetPDGMass();
576 G4double e = mass + 10.0*e 565 G4double e = mass + 10.0*eV;
577 G4double p = std::sqrt(e*e 566 G4double p = std::sqrt(e*e-mass*mass);
578 G4double costheta = 2.0*G4Uniform 567 G4double costheta = 2.0*G4UniformRand() - 1.0;
579 G4double sintheta = std::sqrt((1. 568 G4double sintheta = std::sqrt((1.0 - costheta)*(1.0 + costheta));
580 G4double phi = twopi * G4Uni 569 G4double phi = twopi * G4UniformRand() * rad;
581 G4ThreeVector direction(sintheta*std::cos( 570 G4ThreeVector direction(sintheta*std::cos(phi),sintheta*std::sin(phi),costheta);
582 G4LorentzVector lorentzVector = G4LorentzV 571 G4LorentzVector lorentzVector = G4LorentzVector(direction*p, e);
583 lorentzVector.boost(-boost); 572 lorentzVector.boost(-boost);
584 // Possibility that the following line is not 573 // Possibility that the following line is not correctly carrying over A and Z
585 // from particle definition. Force values. P 574 // from particle definition. Force values. PRT 03/12/2009.
586 // G4Fragment *fragment = 575 // G4Fragment *fragment =
587 // new G4Fragment(lorentzVector, type); 576 // new G4Fragment(lorentzVector, type);
588 G4int A = type->GetBaryonNumber(); 577 G4int A = type->GetBaryonNumber();
589 G4int Z = (G4int) (type->GetPDGCharge() + 578 G4int Z = (G4int) (type->GetPDGCharge() + 1.0E-10);
590 G4Fragment *fragment = 579 G4Fragment *fragment =
591 new G4Fragment(A, Z, lorentzVector); 580 new G4Fragment(A, Z, lorentzVector);
592 if (fragment != nullptr) { fragment->SetCr << 581
593 fragmentVector->push_back(fragment); 582 fragmentVector->push_back(fragment);
594 } 583 }
595 } 584 }
596 ////////////////////////////////////////////// 585 ////////////////////////////////////////////////////////////////////////////////
597 // 586 //
598 void G4WilsonAblationModel::PrintWelcomeMessag 587 void G4WilsonAblationModel::PrintWelcomeMessage ()
599 { 588 {
600 G4cout <<G4endl; 589 G4cout <<G4endl;
601 G4cout <<" ********************************* 590 G4cout <<" *****************************************************************"
602 <<G4endl; 591 <<G4endl;
603 G4cout <<" Nuclear ablation model for nuclea 592 G4cout <<" Nuclear ablation model for nuclear-nuclear interactions activated"
604 <<G4endl; 593 <<G4endl;
605 G4cout <<" (Written by QinetiQ Ltd for the E 594 G4cout <<" (Written by QinetiQ Ltd for the European Space Agency)"
606 <<G4endl; 595 <<G4endl;
607 G4cout <<" !!! WARNING: This model is not we 596 G4cout <<" !!! WARNING: This model is not well validation and should not be used for accurate simulation !!!"
608 <<G4endl; 597 <<G4endl;
609 G4cout <<" ********************************* 598 G4cout <<" *****************************************************************"
610 <<G4endl; 599 <<G4endl;
611 G4cout << G4endl; 600 G4cout << G4endl;
612 601
613 return; 602 return;
614 } 603 }
615 ////////////////////////////////////////////// 604 ////////////////////////////////////////////////////////////////////////////////
616 // 605 //
617 606