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