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

Last change on this file since 1343 was 1340, checked in by garnier, 14 years ago

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