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

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

import all except CVS
File size: 29.1 KB
Line 
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * *
21// * Parts of this code which have been developed by QinetiQ Ltd *
22// * under contract to the European Space Agency (ESA) are the *
23// * intellectual property of ESA. Rights to use, copy, modify and *
24// * redistribute this software for general public use are granted *
25// * in compliance with any licensing, distribution and development *
26// * policy adopted by the Geant4 Collaboration. This code has been *
27// * written by QinetiQ Ltd for the European Space Agency, under ESA *
28// * contract 17191/03/NL/LvH (Aurora Programme). *
29// * *
30// * By using, copying, modifying or distributing the software (or *
31// * any work based on the software) you agree to acknowledge its *
32// * use in resulting scientific publications, and indicate your *
33// * acceptance of all terms of the Geant4 Software license. *
34// ********************************************************************
35//
36// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
37//
38// MODULE: 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//
92G4WilsonAbrasionModel::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//
158G4WilsonAbrasionModel::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//
204G4WilsonAbrasionModel::~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//
215G4HadFinalState *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//
644G4Fragment *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//
748G4double 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//
778void 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//
817void 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//
Note: See TracBrowser for help on using the repository browser.