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