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

Last change on this file since 1326 was 1228, checked in by garnier, 14 years ago

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