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: G4WilsonAbrasionModel.cc |
---|
39 | // |
---|
40 | // Version: B.2 |
---|
41 | // Date: 18/01/05 |
---|
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 | // 6 October 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 | // 18 January 2005, M H Mendenhall, Vanderbilt University, US |
---|
59 | // Pointers to theAbrasionGeometry and products generated by the deexcitation |
---|
60 | // handler deleted to prevent memory leaks. Also particle change of projectile |
---|
61 | // fragment previously not properly defined. |
---|
62 | // |
---|
63 | // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
64 | //////////////////////////////////////////////////////////////////////////////// |
---|
65 | // |
---|
66 | #include "G4WilsonAbrasionModel.hh" |
---|
67 | #include "G4WilsonRadius.hh" |
---|
68 | #include "G4NuclearAbrasionGeometry.hh" |
---|
69 | #include "G4WilsonAblationModel.hh" |
---|
70 | |
---|
71 | #include "G4ExcitationHandler.hh" |
---|
72 | #include "G4Evaporation.hh" |
---|
73 | #include "G4FermiBreakUp.hh" |
---|
74 | #include "G4StatMF.hh" |
---|
75 | #include "G4ParticleDefinition.hh" |
---|
76 | #include "G4DynamicParticle.hh" |
---|
77 | #include "Randomize.hh" |
---|
78 | #include "G4Fragment.hh" |
---|
79 | #include "G4VNuclearDensity.hh" |
---|
80 | #include "G4NuclearShellModelDensity.hh" |
---|
81 | #include "G4NuclearFermiDensity.hh" |
---|
82 | #include "G4FermiMomentum.hh" |
---|
83 | #include "G4ReactionProductVector.hh" |
---|
84 | #include "G4LorentzVector.hh" |
---|
85 | #include "G4ParticleMomentum.hh" |
---|
86 | #include "G4Poisson.hh" |
---|
87 | #include "G4ParticleTable.hh" |
---|
88 | #include "G4IonTable.hh" |
---|
89 | #include "globals.hh" |
---|
90 | //////////////////////////////////////////////////////////////////////////////// |
---|
91 | // |
---|
92 | G4WilsonAbrasionModel::G4WilsonAbrasionModel (G4bool useAblation1) |
---|
93 | :G4HadronicInteraction("G4WilsonAbrasion") |
---|
94 | { |
---|
95 | // |
---|
96 | // |
---|
97 | // Send message to stdout to advise that the G4Abrasion model is being used. |
---|
98 | // |
---|
99 | PrintWelcomeMessage(); |
---|
100 | // |
---|
101 | // |
---|
102 | // Set the default verbose level to 0 - no output. |
---|
103 | // |
---|
104 | verboseLevel = 0; |
---|
105 | useAblation = useAblation1; |
---|
106 | // |
---|
107 | // |
---|
108 | // No de-excitation handler has been supplied - define the default handler. |
---|
109 | // |
---|
110 | theExcitationHandler = new G4ExcitationHandler; |
---|
111 | theExcitationHandlerx = new G4ExcitationHandler; |
---|
112 | if (useAblation) |
---|
113 | { |
---|
114 | theAblation = new G4WilsonAblationModel; |
---|
115 | theAblation->SetVerboseLevel(verboseLevel); |
---|
116 | theExcitationHandler->SetEvaporation(theAblation); |
---|
117 | theExcitationHandlerx->SetEvaporation(theAblation); |
---|
118 | } |
---|
119 | else |
---|
120 | { |
---|
121 | theAblation = NULL; |
---|
122 | G4Evaporation * theEvaporation = new G4Evaporation; |
---|
123 | G4FermiBreakUp * theFermiBreakUp = new G4FermiBreakUp; |
---|
124 | G4StatMF * theMF = new G4StatMF; |
---|
125 | theExcitationHandler->SetEvaporation(theEvaporation); |
---|
126 | theExcitationHandler->SetFermiModel(theFermiBreakUp); |
---|
127 | theExcitationHandler->SetMultiFragmentation(theMF); |
---|
128 | theExcitationHandler->SetMaxAandZForFermiBreakUp(12, 6); |
---|
129 | theExcitationHandler->SetMinEForMultiFrag(5.0*MeV); |
---|
130 | |
---|
131 | theEvaporation = new G4Evaporation; |
---|
132 | theFermiBreakUp = new G4FermiBreakUp; |
---|
133 | theExcitationHandlerx->SetEvaporation(theEvaporation); |
---|
134 | theExcitationHandlerx->SetFermiModel(theFermiBreakUp); |
---|
135 | theExcitationHandlerx->SetMaxAandZForFermiBreakUp(12, 6); |
---|
136 | } |
---|
137 | // |
---|
138 | // |
---|
139 | // Set the minimum and maximum range for the model (despite nomanclature, this |
---|
140 | // is in energy per nucleon number). |
---|
141 | // |
---|
142 | SetMinEnergy(70.0*MeV); |
---|
143 | SetMaxEnergy(10.1*GeV); |
---|
144 | isBlocked = false; |
---|
145 | // |
---|
146 | // |
---|
147 | // npK, when mutiplied by the nuclear Fermi momentum, determines the range of |
---|
148 | // momentum over which the secondary nucleon momentum is sampled. |
---|
149 | // |
---|
150 | npK = 5.0; |
---|
151 | B = 10.0 * MeV; |
---|
152 | third = 1.0 / 3.0; |
---|
153 | conserveEnergy = false; |
---|
154 | conserveMomentum = true; |
---|
155 | } |
---|
156 | //////////////////////////////////////////////////////////////////////////////// |
---|
157 | // |
---|
158 | G4WilsonAbrasionModel::G4WilsonAbrasionModel (G4ExcitationHandler *aExcitationHandler) |
---|
159 | { |
---|
160 | // |
---|
161 | // |
---|
162 | // Send message to stdout to advise that the G4Abrasion model is being used. |
---|
163 | // |
---|
164 | PrintWelcomeMessage(); |
---|
165 | // |
---|
166 | // |
---|
167 | // Set the default verbose level to 0 - no output. |
---|
168 | // |
---|
169 | verboseLevel = 0; |
---|
170 | // |
---|
171 | // |
---|
172 | // The user is able to provide the excitation handler as well as an argument |
---|
173 | // which is provided in this instantiation is used to determine |
---|
174 | // whether the spectators of the interaction are free following the abrasion. |
---|
175 | // |
---|
176 | theExcitationHandler = aExcitationHandler; |
---|
177 | theExcitationHandlerx = new G4ExcitationHandler; |
---|
178 | G4Evaporation * theEvaporation = new G4Evaporation; |
---|
179 | G4FermiBreakUp * theFermiBreakUp = new G4FermiBreakUp; |
---|
180 | theExcitationHandlerx->SetEvaporation(theEvaporation); |
---|
181 | theExcitationHandlerx->SetFermiModel(theFermiBreakUp); |
---|
182 | theExcitationHandlerx->SetMaxAandZForFermiBreakUp(12, 6); |
---|
183 | // |
---|
184 | // |
---|
185 | // Set the minimum and maximum range for the model (despite nomanclature, this |
---|
186 | // is in energy per nucleon number). |
---|
187 | // |
---|
188 | SetMinEnergy(70.0*MeV); |
---|
189 | SetMaxEnergy(10.1*GeV); |
---|
190 | isBlocked = false; |
---|
191 | // |
---|
192 | // |
---|
193 | // npK, when mutiplied by the nuclear Fermi momentum, determines the range of |
---|
194 | // momentum over which the secondary nucleon momentum is sampled. |
---|
195 | // |
---|
196 | npK = 5.0; |
---|
197 | B = 10.0 * MeV; |
---|
198 | third = 1.0 / 3.0; |
---|
199 | conserveEnergy = false; |
---|
200 | conserveMomentum = true; |
---|
201 | } |
---|
202 | //////////////////////////////////////////////////////////////////////////////// |
---|
203 | // |
---|
204 | G4WilsonAbrasionModel::~G4WilsonAbrasionModel () |
---|
205 | { |
---|
206 | // |
---|
207 | // |
---|
208 | // The destructor doesn't have to do a great deal! |
---|
209 | // |
---|
210 | delete theExcitationHandler; |
---|
211 | delete theExcitationHandlerx; |
---|
212 | } |
---|
213 | //////////////////////////////////////////////////////////////////////////////// |
---|
214 | // |
---|
215 | G4HadFinalState *G4WilsonAbrasionModel::ApplyYourself ( |
---|
216 | const G4HadProjectile &theTrack, G4Nucleus &theTarget) |
---|
217 | { |
---|
218 | // |
---|
219 | // |
---|
220 | // The secondaries will be returned in G4HadFinalState &theParticleChange - |
---|
221 | // initialise this. The original track will always be discontinued and |
---|
222 | // secondaries followed. |
---|
223 | // |
---|
224 | theParticleChange.Clear(); |
---|
225 | theParticleChange.SetStatusChange(stopAndKill); |
---|
226 | // |
---|
227 | // |
---|
228 | // Get relevant information about the projectile and target (A, Z, energy/nuc, |
---|
229 | // momentum, etc). |
---|
230 | // |
---|
231 | const G4ParticleDefinition *definitionP = theTrack.GetDefinition(); |
---|
232 | const G4double AP = definitionP->GetBaryonNumber(); |
---|
233 | const G4double ZP = definitionP->GetPDGCharge(); |
---|
234 | G4LorentzVector pP = theTrack.Get4Momentum(); |
---|
235 | G4double E = theTrack.GetKineticEnergy()/AP; |
---|
236 | G4double AT = theTarget.GetN(); |
---|
237 | G4double ZT = theTarget.GetZ(); |
---|
238 | G4double TotalEPre = theTrack.GetTotalEnergy() + |
---|
239 | theTarget.AtomicMass(AT, ZT) + theTarget.GetEnergyDeposit(); |
---|
240 | G4double TotalEPost = 0.0; |
---|
241 | // |
---|
242 | // |
---|
243 | // Determine the radii of the projectile and target nuclei. |
---|
244 | // |
---|
245 | G4WilsonRadius aR; |
---|
246 | G4double rP = aR.GetWilsonRadius(AP); |
---|
247 | G4double rT = aR.GetWilsonRadius(AT); |
---|
248 | G4double rPsq = rP * rP; |
---|
249 | G4double rTsq = rT * rT; |
---|
250 | if (verboseLevel >= 2) |
---|
251 | { |
---|
252 | G4cout <<"########################################" |
---|
253 | <<"########################################" |
---|
254 | <<G4endl; |
---|
255 | G4cout.precision(6); |
---|
256 | G4cout <<"IN G4WilsonAbrasionModel" <<G4endl; |
---|
257 | G4cout <<"Initial projectile A=" <<AP |
---|
258 | <<", Z=" <<ZP |
---|
259 | <<", radius = " <<rP/fermi <<" fm" |
---|
260 | <<G4endl; |
---|
261 | G4cout <<"Initial target A=" <<AT |
---|
262 | <<", Z=" <<ZT |
---|
263 | <<", radius = " <<rT/fermi <<" fm" |
---|
264 | <<G4endl; |
---|
265 | G4cout <<"Projectile momentum and Energy/nuc = " <<pP <<" ," <<E <<G4endl; |
---|
266 | } |
---|
267 | // |
---|
268 | // |
---|
269 | // The following variables are used to determine the impact parameter in the |
---|
270 | // near-field (i.e. taking into consideration the electrostatic repulsion). |
---|
271 | // |
---|
272 | G4double rm = ZP * ZT * elm_coupling / (E * AP); |
---|
273 | G4double r = 0.0; |
---|
274 | G4double rsq = 0.0; |
---|
275 | // |
---|
276 | // |
---|
277 | // Initialise some of the variables which wll be used to calculate the chord- |
---|
278 | // length for nucleons in the projectile and target, and hence calculate the |
---|
279 | // number of abraded nucleons and the excitation energy. |
---|
280 | // |
---|
281 | G4NuclearAbrasionGeometry *theAbrasionGeometry = NULL; |
---|
282 | G4double CT = 0.0; |
---|
283 | G4double F = 0.0; |
---|
284 | G4int Dabr = 0; |
---|
285 | // |
---|
286 | // |
---|
287 | // The following loop is performed until the number of nucleons which are |
---|
288 | // abraded by the process is >1, i.e. an interaction MUST occur. |
---|
289 | // |
---|
290 | while (Dabr == 0) |
---|
291 | { |
---|
292 | // Added by MHM 20050119 to fix leaking memory on second pass through this loop |
---|
293 | if (theAbrasionGeometry) |
---|
294 | { |
---|
295 | delete theAbrasionGeometry; |
---|
296 | theAbrasionGeometry = NULL; |
---|
297 | } |
---|
298 | // |
---|
299 | // |
---|
300 | // Sample the impact parameter. For the moment, this class takes account of |
---|
301 | // electrostatic effects on the impact parameter, but (like HZETRN AND NUCFRG2) |
---|
302 | // does not make any correction for the effects of nuclear-nuclear repulsion. |
---|
303 | // |
---|
304 | G4double rPT = rP + rT; |
---|
305 | G4double rPTsq = rPT * rPT; |
---|
306 | r = 1.1 * rPT; |
---|
307 | while (r > rPT) |
---|
308 | { |
---|
309 | G4double bsq = rPTsq * G4UniformRand(); |
---|
310 | r = (rm + std::sqrt(rm*rm + 4.0*bsq)) / 2.0; |
---|
311 | } |
---|
312 | rsq = r * r; |
---|
313 | // |
---|
314 | // |
---|
315 | // Now determine the chord-length through the target nucleus. |
---|
316 | // |
---|
317 | if (rT > rP) |
---|
318 | { |
---|
319 | G4double x = (rPsq + rsq - rTsq) / 2.0 / r; |
---|
320 | if (x > 0.0) CT = 2.0 * std::sqrt(rTsq - x*x); |
---|
321 | else CT = 2.0 * std::sqrt(rTsq - rsq); |
---|
322 | } |
---|
323 | else |
---|
324 | { |
---|
325 | G4double x = (rTsq + rsq - rPsq) / 2.0 / r; |
---|
326 | if (x > 0.0) CT = 2.0 * std::sqrt(rTsq - x*x); |
---|
327 | else CT = 2.0 * rT; |
---|
328 | } |
---|
329 | // |
---|
330 | // |
---|
331 | // Determine the number of abraded nucleons. Note that the mean number of |
---|
332 | // abraded nucleons is used to sample the Poisson distribution. The Poisson |
---|
333 | // distribution is sampled only ten times with the current impact parameter, |
---|
334 | // and if it fails after this to find a case for which the number of abraded |
---|
335 | // nucleons >1, the impact parameter is re-sampled. |
---|
336 | // |
---|
337 | theAbrasionGeometry = new G4NuclearAbrasionGeometry(AP,AT,r); |
---|
338 | F = theAbrasionGeometry->F(); |
---|
339 | G4double lambda = 16.6*fermi / std::pow(E/MeV,0.26); |
---|
340 | G4double Mabr = F * AP * (1.0 - std::exp(-CT/lambda)); |
---|
341 | G4long n = 0; |
---|
342 | for (G4int i = 0; i<10; i++) |
---|
343 | { |
---|
344 | n = G4Poisson(Mabr); |
---|
345 | if (n > 0) |
---|
346 | { |
---|
347 | if (n>AP) Dabr = (G4int) AP; |
---|
348 | else Dabr = (G4int) n; |
---|
349 | break; |
---|
350 | } |
---|
351 | } |
---|
352 | } |
---|
353 | if (verboseLevel >= 2) |
---|
354 | { |
---|
355 | G4cout <<G4endl; |
---|
356 | G4cout <<"Impact parameter = " <<r/fermi <<" fm" <<G4endl; |
---|
357 | G4cout <<"# Abraded nucleons = " <<Dabr <<G4endl; |
---|
358 | } |
---|
359 | // |
---|
360 | // |
---|
361 | // The number of abraded nucleons must be no greater than the number of |
---|
362 | // nucleons in either the projectile or the target. If AP - Dabr < 2 or |
---|
363 | // AT - Dabr < 2 then either we have only a nucleon left behind in the |
---|
364 | // projectile/target or we've tried to abrade too many nucleons - and Dabr |
---|
365 | // should be limited. |
---|
366 | // |
---|
367 | if (AP - (G4double) Dabr < 2.0) Dabr = (G4int) AP; |
---|
368 | if (AT - (G4double) Dabr < 2.0) Dabr = (G4int) AT; |
---|
369 | // |
---|
370 | // |
---|
371 | // Determine the abraded secondary nucleons from the projectile. *fragmentP |
---|
372 | // is a pointer to the prefragment from the projectile and nSecP is the number |
---|
373 | // of nucleons in theParticleChange which have been abraded. The total energy |
---|
374 | // from these is determined. |
---|
375 | // |
---|
376 | G4ThreeVector boost = pP.findBoostToCM(); |
---|
377 | G4Fragment *fragmentP = GetAbradedNucleons (Dabr, AP, ZP, rP); |
---|
378 | G4int nSecP = theParticleChange.GetNumberOfSecondaries(); |
---|
379 | G4int i = 0; |
---|
380 | for (i=0; i<nSecP; i++) |
---|
381 | { |
---|
382 | TotalEPost += theParticleChange.GetSecondary(i)-> |
---|
383 | GetParticle()->GetTotalEnergy(); |
---|
384 | } |
---|
385 | // |
---|
386 | // |
---|
387 | // Determine the number of spectators in the interaction region for the |
---|
388 | // projectile. |
---|
389 | // |
---|
390 | G4int DspcP = (G4int) (AP*F) - Dabr; |
---|
391 | if (DspcP <= 0) DspcP = 0; |
---|
392 | else if (DspcP > AP-Dabr) DspcP = ((G4int) AP) - Dabr; |
---|
393 | // |
---|
394 | // |
---|
395 | // Determine excitation energy associated with excess surface area of the |
---|
396 | // projectile (EsP) and the excitation due to scattering of nucleons which are |
---|
397 | // retained within the projectile (ExP). Add the total energy from the excited |
---|
398 | // nucleus to the total energy of the secondaries. |
---|
399 | // |
---|
400 | G4bool excitationAbsorbedByProjectile = false; |
---|
401 | if (fragmentP != NULL) |
---|
402 | { |
---|
403 | G4double EsP = theAbrasionGeometry->GetExcitationEnergyOfProjectile(); |
---|
404 | G4double ExP = 0.0; |
---|
405 | if (Dabr < AT) |
---|
406 | excitationAbsorbedByProjectile = G4UniformRand() < 0.5; |
---|
407 | if (excitationAbsorbedByProjectile) |
---|
408 | ExP = GetNucleonInducedExcitation(rP, rT, r); |
---|
409 | G4double xP = EsP + ExP; |
---|
410 | if (xP > B*(AP-Dabr)) xP = B*(AP-Dabr); |
---|
411 | G4LorentzVector lorentzVector = fragmentP->GetMomentum(); |
---|
412 | lorentzVector.setE(lorentzVector.e()+xP); |
---|
413 | fragmentP->SetMomentum(lorentzVector); |
---|
414 | TotalEPost += lorentzVector.e(); |
---|
415 | } |
---|
416 | G4double EMassP = TotalEPost; |
---|
417 | // |
---|
418 | // |
---|
419 | // Determine the abraded secondary nucleons from the target. Note that it's |
---|
420 | // assumed that the same number of nucleons are abraded from the target as for |
---|
421 | // the projectile, and obviously no boost is applied to the products. *fragmentT |
---|
422 | // is a pointer to the prefragment from the target and nSec is the total number |
---|
423 | // of nucleons in theParticleChange which have been abraded. The total energy |
---|
424 | // from these is determined. |
---|
425 | // |
---|
426 | G4Fragment *fragmentT = GetAbradedNucleons (Dabr, AT, ZT, rT); |
---|
427 | G4int nSec = theParticleChange.GetNumberOfSecondaries(); |
---|
428 | for (i=nSecP; i<nSec; i++) |
---|
429 | { |
---|
430 | TotalEPost += theParticleChange.GetSecondary(i)-> |
---|
431 | GetParticle()->GetTotalEnergy(); |
---|
432 | } |
---|
433 | // |
---|
434 | // |
---|
435 | // Determine the number of spectators in the interaction region for the |
---|
436 | // target. |
---|
437 | // |
---|
438 | G4int DspcT = (G4int) (AT*F) - Dabr; |
---|
439 | if (DspcT <= 0) DspcT = 0; |
---|
440 | else if (DspcT > AP-Dabr) DspcT = ((G4int) AT) - Dabr; |
---|
441 | // |
---|
442 | // |
---|
443 | // Determine excitation energy associated with excess surface area of the |
---|
444 | // target (EsT) and the excitation due to scattering of nucleons which are |
---|
445 | // retained within the target (ExT). Add the total energy from the excited |
---|
446 | // nucleus to the total energy of the secondaries. |
---|
447 | // |
---|
448 | if (fragmentT != NULL) |
---|
449 | { |
---|
450 | G4double EsT = theAbrasionGeometry->GetExcitationEnergyOfTarget(); |
---|
451 | G4double ExT = 0.0; |
---|
452 | if (!excitationAbsorbedByProjectile) |
---|
453 | ExT = GetNucleonInducedExcitation(rT, rP, r); |
---|
454 | G4double xT = EsT + ExT; |
---|
455 | if (xT > B*(AT-Dabr)) xT = B*(AT-Dabr); |
---|
456 | G4LorentzVector lorentzVector = fragmentT->GetMomentum(); |
---|
457 | lorentzVector.setE(lorentzVector.e()+xT); |
---|
458 | fragmentT->SetMomentum(lorentzVector); |
---|
459 | TotalEPost += lorentzVector.e(); |
---|
460 | } |
---|
461 | // |
---|
462 | // |
---|
463 | // Now determine the difference between the pre and post interaction |
---|
464 | // energy - this will be used to determine the Lorentz boost if conservation |
---|
465 | // of energy is to be imposed/attempted. |
---|
466 | // |
---|
467 | G4double deltaE = TotalEPre - TotalEPost; |
---|
468 | if (deltaE > 0.0 && conserveEnergy) |
---|
469 | { |
---|
470 | G4double beta = std::sqrt(1.0 - EMassP*EMassP/std::pow(deltaE+EMassP,2.0)); |
---|
471 | boost = boost / boost.mag() * beta; |
---|
472 | } |
---|
473 | // |
---|
474 | // |
---|
475 | // Now boost the secondaries from the projectile. |
---|
476 | // |
---|
477 | G4ThreeVector pBalance = pP.vect(); |
---|
478 | for (i=0; i<nSecP; i++) |
---|
479 | { |
---|
480 | G4DynamicParticle *dynamicP = theParticleChange.GetSecondary(i)-> |
---|
481 | GetParticle(); |
---|
482 | G4LorentzVector lorentzVector = dynamicP->Get4Momentum(); |
---|
483 | lorentzVector.boost(-boost); |
---|
484 | dynamicP->Set4Momentum(lorentzVector); |
---|
485 | pBalance -= lorentzVector.vect(); |
---|
486 | } |
---|
487 | // |
---|
488 | // |
---|
489 | // Set the boost for the projectile prefragment. This is now based on the |
---|
490 | // conservation of momentum. However, if the user selected momentum of the |
---|
491 | // prefragment is not to be conserved this simply boosted to the velocity of the |
---|
492 | // original projectile times the ratio of the unexcited to the excited mass |
---|
493 | // of the prefragment (the excitation increases the effective mass of the |
---|
494 | // prefragment, and therefore modifying the boost is an attempt to prevent |
---|
495 | // the momentum of the prefragment being excessive). |
---|
496 | // |
---|
497 | if (fragmentP != NULL) |
---|
498 | { |
---|
499 | G4LorentzVector lorentzVector = fragmentP->GetMomentum(); |
---|
500 | G4double m = lorentzVector.m(); |
---|
501 | if (conserveMomentum) |
---|
502 | fragmentP->SetMomentum |
---|
503 | (G4LorentzVector(pBalance,std::sqrt(pBalance.mag2()+m*m+1.0*eV*eV))); |
---|
504 | else |
---|
505 | { |
---|
506 | G4double mg = fragmentP->GetGroundStateMass(); |
---|
507 | fragmentP->SetMomentum(lorentzVector.boost(-boost * mg/m)); |
---|
508 | } |
---|
509 | } |
---|
510 | // |
---|
511 | // |
---|
512 | // Output information to user if verbose information requested. |
---|
513 | // |
---|
514 | if (verboseLevel >= 2) |
---|
515 | { |
---|
516 | G4cout <<G4endl; |
---|
517 | G4cout <<"-----------------------------------" <<G4endl; |
---|
518 | G4cout <<"Secondary nucleons from projectile:" <<G4endl; |
---|
519 | G4cout <<"-----------------------------------" <<G4endl; |
---|
520 | G4cout.precision(7); |
---|
521 | for (i=0; i<nSecP; i++) |
---|
522 | { |
---|
523 | G4cout <<"Particle # " <<i <<G4endl; |
---|
524 | theParticleChange.GetSecondary(i)->GetParticle()->DumpInfo(); |
---|
525 | G4DynamicParticle *dyn = theParticleChange.GetSecondary(i)->GetParticle(); |
---|
526 | G4cout <<"New nucleon (P) " <<dyn->GetDefinition()->GetParticleName() |
---|
527 | <<" : " <<dyn->Get4Momentum() |
---|
528 | <<G4endl; |
---|
529 | } |
---|
530 | G4cout <<"---------------------------" <<G4endl; |
---|
531 | G4cout <<"The projectile prefragment:" <<G4endl; |
---|
532 | G4cout <<"---------------------------" <<G4endl; |
---|
533 | if (fragmentP != NULL) |
---|
534 | G4cout <<*fragmentP <<G4endl; |
---|
535 | else |
---|
536 | G4cout <<"(No residual prefragment)" <<G4endl; |
---|
537 | G4cout <<G4endl; |
---|
538 | G4cout <<"-------------------------------" <<G4endl; |
---|
539 | G4cout <<"Secondary nucleons from target:" <<G4endl; |
---|
540 | G4cout <<"-------------------------------" <<G4endl; |
---|
541 | G4cout.precision(7); |
---|
542 | for (i=nSecP; i<nSec; i++) |
---|
543 | { |
---|
544 | G4cout <<"Particle # " <<i <<G4endl; |
---|
545 | theParticleChange.GetSecondary(i)->GetParticle()->DumpInfo(); |
---|
546 | G4DynamicParticle *dyn = theParticleChange.GetSecondary(i)->GetParticle(); |
---|
547 | G4cout <<"New nucleon (T) " <<dyn->GetDefinition()->GetParticleName() |
---|
548 | <<" : " <<dyn->Get4Momentum() |
---|
549 | <<G4endl; |
---|
550 | } |
---|
551 | G4cout <<"-----------------------" <<G4endl; |
---|
552 | G4cout <<"The target prefragment:" <<G4endl; |
---|
553 | G4cout <<"-----------------------" <<G4endl; |
---|
554 | if (fragmentT != NULL) |
---|
555 | G4cout <<*fragmentT <<G4endl; |
---|
556 | else |
---|
557 | G4cout <<"(No residual prefragment)" <<G4endl; |
---|
558 | } |
---|
559 | // |
---|
560 | // |
---|
561 | // Now we can decay the nuclear fragments if present. The secondaries are |
---|
562 | // collected and boosted as well. This is performed first for the projectile... |
---|
563 | // |
---|
564 | if (fragmentP !=NULL) |
---|
565 | { |
---|
566 | G4ReactionProductVector *products = NULL; |
---|
567 | if (fragmentP->GetZ() != fragmentP->GetA()) |
---|
568 | products = theExcitationHandler->BreakItUp(*fragmentP); |
---|
569 | else |
---|
570 | products = theExcitationHandlerx->BreakItUp(*fragmentP); |
---|
571 | delete fragmentP; |
---|
572 | fragmentP = NULL; |
---|
573 | |
---|
574 | G4ReactionProductVector::iterator iter; |
---|
575 | for (iter = products->begin(); iter != products->end(); ++iter) |
---|
576 | { |
---|
577 | G4DynamicParticle *secondary = |
---|
578 | new G4DynamicParticle((*iter)->GetDefinition(), |
---|
579 | (*iter)->GetTotalEnergy(), (*iter)->GetMomentum()); |
---|
580 | theParticleChange.AddSecondary (secondary); // Added MHM 20050118 |
---|
581 | G4String particleName = (*iter)->GetDefinition()->GetParticleName(); |
---|
582 | delete (*iter); // get rid of leftover particle def! // Added MHM 20050118 |
---|
583 | if (verboseLevel >= 2 && particleName.find("[",0) < particleName.size()) |
---|
584 | { |
---|
585 | G4cout <<"------------------------" <<G4endl; |
---|
586 | G4cout <<"The projectile fragment:" <<G4endl; |
---|
587 | G4cout <<"------------------------" <<G4endl; |
---|
588 | G4cout <<" fragmentP = " <<particleName |
---|
589 | <<" Energy = " <<secondary->GetKineticEnergy() |
---|
590 | <<G4endl; |
---|
591 | } |
---|
592 | } |
---|
593 | delete products; // Added MHM 20050118 |
---|
594 | } |
---|
595 | // |
---|
596 | // |
---|
597 | // Now decay the target nucleus - no boost is applied since in this |
---|
598 | // approximation it is assumed that there is negligible momentum transfer from |
---|
599 | // the projectile. |
---|
600 | // |
---|
601 | if (fragmentT != NULL) |
---|
602 | { |
---|
603 | G4ReactionProductVector *products = NULL; |
---|
604 | if (fragmentT->GetZ() != fragmentT->GetA()) |
---|
605 | products = theExcitationHandler->BreakItUp(*fragmentT); |
---|
606 | else |
---|
607 | products = theExcitationHandlerx->BreakItUp(*fragmentT); |
---|
608 | delete fragmentT; |
---|
609 | fragmentT = NULL; |
---|
610 | |
---|
611 | G4ReactionProductVector::iterator iter; |
---|
612 | for (iter = products->begin(); iter != products->end(); ++iter) |
---|
613 | { |
---|
614 | G4DynamicParticle *secondary = |
---|
615 | new G4DynamicParticle((*iter)->GetDefinition(), |
---|
616 | (*iter)->GetTotalEnergy(), (*iter)->GetMomentum()); |
---|
617 | theParticleChange.AddSecondary (secondary); |
---|
618 | G4String particleName = (*iter)->GetDefinition()->GetParticleName(); |
---|
619 | delete (*iter); // get rid of leftover particle def! // Added MHM 20050118 |
---|
620 | if (verboseLevel >= 2 && particleName.find("[",0) < particleName.size()) |
---|
621 | { |
---|
622 | G4cout <<"--------------------" <<G4endl; |
---|
623 | G4cout <<"The target fragment:" <<G4endl; |
---|
624 | G4cout <<"--------------------" <<G4endl; |
---|
625 | G4cout <<" fragmentT = " <<particleName |
---|
626 | <<" Energy = " <<secondary->GetKineticEnergy() |
---|
627 | <<G4endl; |
---|
628 | } |
---|
629 | } |
---|
630 | delete products; // Added MHM 20050118 |
---|
631 | } |
---|
632 | |
---|
633 | if (verboseLevel >= 2) |
---|
634 | G4cout <<"########################################" |
---|
635 | <<"########################################" |
---|
636 | <<G4endl; |
---|
637 | |
---|
638 | delete theAbrasionGeometry; |
---|
639 | |
---|
640 | return &theParticleChange; |
---|
641 | } |
---|
642 | //////////////////////////////////////////////////////////////////////////////// |
---|
643 | // |
---|
644 | G4Fragment *G4WilsonAbrasionModel::GetAbradedNucleons (G4int Dabr, G4double A, |
---|
645 | G4double Z, G4double r) |
---|
646 | { |
---|
647 | // |
---|
648 | // |
---|
649 | // Initialise variables. tau is the Fermi radius of the nucleus. The variables |
---|
650 | // p..., C... and g(amma) are used to help sample the secondary nucleon |
---|
651 | // spectrum. |
---|
652 | // |
---|
653 | |
---|
654 | G4double pK = hbarc * std::pow(9.0 * pi / 4.0 * A, third) / (1.29 * r); |
---|
655 | if (A <= 24.0) pK *= -0.229*std::pow(A,third) + 1.62; |
---|
656 | G4double pKsq = pK * pK; |
---|
657 | G4double p1sq = 2.0/5.0 * pKsq; |
---|
658 | G4double p2sq = 6.0/5.0 * pKsq; |
---|
659 | G4double p3sq = 500.0 * 500.0; |
---|
660 | G4double C1 = 1.0; |
---|
661 | G4double C2 = 0.03; |
---|
662 | G4double C3 = 0.0002; |
---|
663 | G4double g = 90.0 * MeV; |
---|
664 | G4double maxn = C1 + C2 + C3; |
---|
665 | // |
---|
666 | // |
---|
667 | // initialise the number of secondary nucleons abraded to zero, and initially set |
---|
668 | // the type of nucleon abraded to proton ... just for now. |
---|
669 | // |
---|
670 | G4double Aabr = 0.0; |
---|
671 | G4double Zabr = 0.0; |
---|
672 | G4ParticleDefinition *typeNucleon = G4Proton::ProtonDefinition(); |
---|
673 | G4DynamicParticle *dynamicNucleon = NULL; |
---|
674 | G4ParticleMomentum pabr(0.0, 0.0, 0.0); |
---|
675 | // |
---|
676 | // |
---|
677 | // Now go through each abraded nucleon and sample type, spectrum and angle. |
---|
678 | // |
---|
679 | for (G4int i=0; i<Dabr; i++) |
---|
680 | { |
---|
681 | // |
---|
682 | // |
---|
683 | // Sample the nucleon momentum distribution by simple rejection techniques. We |
---|
684 | // reject values of p == 0.0 since this causes bad behaviour in the sinh term. |
---|
685 | // |
---|
686 | G4double p = 0.0; |
---|
687 | G4bool found = false; |
---|
688 | while (!found) |
---|
689 | { |
---|
690 | while (p <= 0.0) p = npK * pK * G4UniformRand(); |
---|
691 | G4double psq = p * p; |
---|
692 | found = maxn * G4UniformRand() < C1*std::exp(-psq/p1sq/2.0) + |
---|
693 | C2*std::exp(-psq/p2sq/2.0) + C3*std::exp(-psq/p3sq/2.0) + p/g/std::sinh(p/g); |
---|
694 | } |
---|
695 | // |
---|
696 | // |
---|
697 | // Determine the type of particle abraded. Can only be proton or neutron, |
---|
698 | // and the probability is determine to be proportional to the ratio as found |
---|
699 | // in the nucleus at each stage. |
---|
700 | // |
---|
701 | G4double prob = (Z-Zabr)/(A-Aabr); |
---|
702 | if (G4UniformRand()<prob) |
---|
703 | { |
---|
704 | Zabr++; |
---|
705 | typeNucleon = G4Proton::ProtonDefinition(); |
---|
706 | } |
---|
707 | else |
---|
708 | typeNucleon = G4Neutron::NeutronDefinition(); |
---|
709 | Aabr++; |
---|
710 | // |
---|
711 | // |
---|
712 | // The angular distribution of the secondary nucleons is approximated to an |
---|
713 | // isotropic distribution in the rest frame of the nucleus (this will be Lorentz |
---|
714 | // boosted later. |
---|
715 | // |
---|
716 | G4double costheta = 2.*G4UniformRand()-1.0; |
---|
717 | G4double sintheta = std::sqrt((1.0 - costheta)*(1.0 + costheta)); |
---|
718 | G4double phi = 2.0*pi*G4UniformRand()*rad; |
---|
719 | G4ThreeVector direction(sintheta*std::cos(phi),sintheta*std::sin(phi),costheta); |
---|
720 | G4double nucleonMass = typeNucleon->GetPDGMass(); |
---|
721 | G4double E = std::sqrt(p*p + nucleonMass*nucleonMass)-nucleonMass; |
---|
722 | dynamicNucleon = new G4DynamicParticle(typeNucleon,direction,E); |
---|
723 | theParticleChange.AddSecondary (dynamicNucleon); |
---|
724 | pabr += p*direction; |
---|
725 | } |
---|
726 | // |
---|
727 | // |
---|
728 | // Next determine the details of the nuclear prefragment .. that is if there |
---|
729 | // is one or more protons in the residue. (Note that the 1 eV in the total |
---|
730 | // energy is a safety factor to avoid any possibility of negative rest mass |
---|
731 | // energy.) |
---|
732 | // |
---|
733 | G4Fragment *fragment = NULL; |
---|
734 | if (Z-Zabr>=1.0) |
---|
735 | { |
---|
736 | G4double ionMass = G4ParticleTable::GetParticleTable()->GetIonTable()-> |
---|
737 | GetIonMass(G4lrint(Z-Zabr),G4lrint(A-Aabr)); |
---|
738 | G4double E = std::sqrt(pabr.mag2() + ionMass*ionMass); |
---|
739 | G4LorentzVector lorentzVector = G4LorentzVector(-pabr, E + 1.0*eV); |
---|
740 | fragment = |
---|
741 | new G4Fragment((G4int) (A-Aabr), (G4int) (Z-Zabr), lorentzVector); |
---|
742 | } |
---|
743 | |
---|
744 | return fragment; |
---|
745 | } |
---|
746 | //////////////////////////////////////////////////////////////////////////////// |
---|
747 | // |
---|
748 | G4double G4WilsonAbrasionModel::GetNucleonInducedExcitation |
---|
749 | (G4double rP, G4double rT, G4double r) |
---|
750 | { |
---|
751 | // |
---|
752 | // |
---|
753 | // Initialise variables. |
---|
754 | // |
---|
755 | G4double Cl = 0.0; |
---|
756 | G4double rPsq = rP * rP; |
---|
757 | G4double rTsq = rT * rT; |
---|
758 | G4double rsq = r * r; |
---|
759 | // |
---|
760 | // |
---|
761 | // Depending upon the impact parameter, a different form of the chord length is |
---|
762 | // is used. |
---|
763 | // |
---|
764 | if (r > rT) Cl = 2.0*std::sqrt(rPsq + 2.0*r*rT - rsq - rTsq); |
---|
765 | else Cl = 2.0*rP; |
---|
766 | |
---|
767 | G4double bP = (rPsq+rsq-rTsq)/2.0/r; |
---|
768 | G4double Ct = 2.0*std::sqrt(rPsq - bP*bP); |
---|
769 | |
---|
770 | G4double Ex = 13.0 * Cl / fermi; |
---|
771 | if (Ct > 1.5*fermi) |
---|
772 | Ex += 13.0 * Cl / fermi /3.0 * (Ct/fermi - 1.5); |
---|
773 | |
---|
774 | return Ex; |
---|
775 | } |
---|
776 | //////////////////////////////////////////////////////////////////////////////// |
---|
777 | // |
---|
778 | void G4WilsonAbrasionModel::SetUseAblation (G4bool useAblation1) |
---|
779 | { |
---|
780 | if (useAblation != useAblation1) |
---|
781 | { |
---|
782 | useAblation = useAblation1; |
---|
783 | delete theExcitationHandler; |
---|
784 | delete theExcitationHandlerx; |
---|
785 | theExcitationHandler = new G4ExcitationHandler; |
---|
786 | theExcitationHandlerx = new G4ExcitationHandler; |
---|
787 | if (useAblation) |
---|
788 | { |
---|
789 | theAblation = new G4WilsonAblationModel; |
---|
790 | theAblation->SetVerboseLevel(verboseLevel); |
---|
791 | theExcitationHandler->SetEvaporation(theAblation); |
---|
792 | theExcitationHandlerx->SetEvaporation(theAblation); |
---|
793 | } |
---|
794 | else |
---|
795 | { |
---|
796 | theAblation = NULL; |
---|
797 | G4Evaporation * theEvaporation = new G4Evaporation; |
---|
798 | G4FermiBreakUp * theFermiBreakUp = new G4FermiBreakUp; |
---|
799 | G4StatMF * theMF = new G4StatMF; |
---|
800 | theExcitationHandler->SetEvaporation(theEvaporation); |
---|
801 | theExcitationHandler->SetFermiModel(theFermiBreakUp); |
---|
802 | theExcitationHandler->SetMultiFragmentation(theMF); |
---|
803 | theExcitationHandler->SetMaxAandZForFermiBreakUp(12, 6); |
---|
804 | theExcitationHandler->SetMinEForMultiFrag(5.0*MeV); |
---|
805 | |
---|
806 | theEvaporation = new G4Evaporation; |
---|
807 | theFermiBreakUp = new G4FermiBreakUp; |
---|
808 | theExcitationHandlerx->SetEvaporation(theEvaporation); |
---|
809 | theExcitationHandlerx->SetFermiModel(theFermiBreakUp); |
---|
810 | theExcitationHandlerx->SetMaxAandZForFermiBreakUp(12, 6); |
---|
811 | } |
---|
812 | } |
---|
813 | return; |
---|
814 | } |
---|
815 | //////////////////////////////////////////////////////////////////////////////// |
---|
816 | // |
---|
817 | void G4WilsonAbrasionModel::PrintWelcomeMessage () |
---|
818 | { |
---|
819 | G4cout <<G4endl; |
---|
820 | G4cout <<" *****************************************************************" |
---|
821 | <<G4endl; |
---|
822 | G4cout <<" Nuclear abrasion model for nuclear-nuclear interactions activated" |
---|
823 | <<G4endl; |
---|
824 | G4cout <<" (Written by QinetiQ Ltd for the European Space Agency)" |
---|
825 | <<G4endl; |
---|
826 | G4cout <<" *****************************************************************" |
---|
827 | <<G4endl; |
---|
828 | G4cout << G4endl; |
---|
829 | |
---|
830 | return; |
---|
831 | } |
---|
832 | //////////////////////////////////////////////////////////////////////////////// |
---|
833 | // |
---|