| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // *
21 // * Parts of this code which have been devel
22 // * under contract to the European Space Agen
23 // * intellectual property of ESA. Rights to u
24 // * redistribute this software for general pu
25 // * in compliance with any licensing, distrib
26 // * policy adopted by the Geant4 Collaboratio
27 // * written by QinetiQ Ltd for the European S
28 // * contract 17191/03/NL/LvH (Aurora Programm
29 // *
30 // * By using, copying, modifying or distri
31 // * any work based on the software) you ag
32 // * use in resulting scientific publicati
33 // * acceptance of all terms of the Geant4 Sof
34 // *******************************************
35 //
36 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
37 //
38 // MODULE: G4WilsonAbrasionModel.
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,
53 // Created.
54 //
55 // 15 March 2004, P R Truscott, QinetiQ Ltd, U
56 // Beta release
57 //
58 // 18 January 2005, M H Mendenhall, Vanderbilt
59 // Pointers to theAbrasionGeometry and product
60 // handler deleted to prevent memory leaks. A
61 // fragment previously not properly defined.
62 //
63 // 08 December 2009, P R Truscott, QinetiQ Ltd
64 // ver 1.0
65 // There was originally a possibility of the m
66 // considering Couloumb repulsion, rm, being t
67 // rm >= fradius * (rP + rT)
68 // where fradius is currently 0.99, then it is
69 // unchanged. Additional conditions to escape
70 // parameter: if the loop counter evtcnt reach
71 // continued.
72 // Additional clauses have been included in
73 // G4WilsonAbrasionModel::GetNucleonInduced
74 // Previously it was possible to get sqrt of n
75 // algorithm not properly defined if either:
76 // rT > rP && rsq < rTsq - rPsq) or (rP > r
77 //
78 // 12 June 2012, A. Ribon, CERN, Switzerland
79 // Fixing trivial warning errors of shadowed v
80 //
81 // 4 August 2015, A. Ribon, CERN, Switzerland
82 // Replacing std::exp and std::pow with the fa
83 //
84 // 7 August 2015, A. Ribon, CERN, Switzerland
85 // Checking of 'while' loops.
86 //
87 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
88 //////////////////////////////////////////////
89
90 #include "G4WilsonAbrasionModel.hh"
91 #include "G4WilsonRadius.hh"
92 #include "G4NuclearAbrasionGeometry.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"
103 #include "G4ReactionProductVector.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 #include "G4PhysicsModelCatalog.hh"
115
116
117 G4WilsonAbrasionModel::G4WilsonAbrasionModel(G
118 : G4HadronicInteraction("G4WilsonAbrasion"),
119 {
120 // Send message to stdout to advise that the
121 PrintWelcomeMessage();
122
123 // Set the default verbose level to 0 - no o
124 verboseLevel = 0;
125 useAblation = useAblation1;
126 theAblation = nullptr;
127
128 // No de-excitation handler has been supplie
129
130 theExcitationHandler = new G4ExcitationHandl
131 if (useAblation)
132 {
133 theAblation = new G4WilsonAblationModel;
134 theAblation->SetVerboseLevel(verboseLevel)
135 theExcitationHandler->SetEvaporation(theAb
136 }
137
138 // Set the minimum and maximum range for the
139 // this is in energy per nucleon number).
140
141 SetMinEnergy(70.0*MeV);
142 SetMaxEnergy(10.1*GeV);
143 isBlocked = false;
144
145 // npK, when mutiplied by the nuclear Fermi
146 // momentum over which the secondary nucleon
147
148 r0sq = 0.0;
149 npK = 5.0;
150 B = 10.0 * MeV;
151 third = 1.0 / 3.0;
152 fradius = 0.99;
153 conserveEnergy = false;
154 conserveMomentum = true;
155
156 // Creator model ID for the secondaries crea
157 secID = G4PhysicsModelCatalog::GetModelID( "
158 }
159
160 void G4WilsonAbrasionModel::ModelDescription(s
161 {
162 outFile << "G4WilsonAbrasionModel is a macro
163 << "nucleus-nucleus collisions using
164 << "The smaller projectile nucleus g
165 << "target nucleus, leaving a residu
166 << "region where the projectile and
167 << "is then treated as a highly exci
168 << "model is based on the NUCFRG2 mo
169 << "projectile energies between 70 M
170 }
171
172 G4WilsonAbrasionModel::G4WilsonAbrasionModel(G
173 G4HadronicInteraction("G4WilsonAbrasion"), se
174 {
175 // Send message to stdout to advise that the G
176
177 PrintWelcomeMessage();
178
179 // Set the default verbose level to 0 - no out
180
181 verboseLevel = 0;
182
183 theAblation = nullptr; //A.R. 26-Jul-2012
184 useAblation = false; //A.R. 14-Aug-2012
185
186 //
187 // The user is able to provide the excitation
188 // which is provided in this instantiation is
189 // whether the spectators of the interaction a
190 //
191 theExcitationHandler = aExcitationHandler;
192 //
193 //
194 // Set the minimum and maximum range for the m
195 // is in energy per nucleon number).
196 //
197 SetMinEnergy(70.0*MeV);
198 SetMaxEnergy(10.1*GeV);
199 isBlocked = false;
200 //
201 //
202 // npK, when mutiplied by the nuclear Fermi mo
203 // momentum over which the secondary nucleon m
204 //
205 r0sq = 0.0;
206 npK = 5.0;
207 B = 10.0 * MeV;
208 third = 1.0 / 3.0;
209 fradius = 0.99;
210 conserveEnergy = false;
211 conserveMomentum = true;
212
213 // Creator model ID for the secondaries crea
214 secID = G4PhysicsModelCatalog::GetModelID( "
215 }
216 //////////////////////////////////////////////
217 //
218 G4WilsonAbrasionModel::~G4WilsonAbrasionModel(
219 {
220 delete theExcitationHandler;
221 }
222 //////////////////////////////////////////////
223 //
224 G4HadFinalState *G4WilsonAbrasionModel::ApplyY
225 const G4HadProjectile &theTrack, G4Nucleus &
226 {
227 //
228 //
229 // The secondaries will be returned in G4HadFi
230 // initialise this. The original track will a
231 // secondaries followed.
232 //
233 theParticleChange.Clear();
234 theParticleChange.SetStatusChange(stopAndKil
235 //
236 //
237 // Get relevant information about the projecti
238 // momentum, etc).
239 //
240 const G4ParticleDefinition *definitionP = th
241 const G4double AP = definitionP->GetBaryonN
242 const G4double ZP = definitionP->GetPDGChar
243 G4LorentzVector pP = theTrack.Get4Momentum()
244 G4double E = theTrack.GetKineticEner
245 G4double AT = theTarget.GetA_asInt();
246 G4double ZT = theTarget.GetZ_asInt();
247 G4double TotalEPre = theTrack.GetTotalEnergy
248 theTarget.AtomicMass(AT, ZT) + theTarget.G
249 G4double TotalEPost = 0.0;
250 //
251 //
252 // Determine the radii of the projectile and t
253 //
254 G4WilsonRadius aR;
255 G4double rP = aR.GetWilsonRadius(AP);
256 G4double rT = aR.GetWilsonRadius(AT);
257 G4double rPsq = rP * rP;
258 G4double rTsq = rT * rT;
259 if (verboseLevel >= 2)
260 {
261 G4cout <<"################################
262 <<"################################
263 <<G4endl;
264 G4cout.precision(6);
265 G4cout <<"IN G4WilsonAbrasionModel" <<G4en
266 G4cout <<"Initial projectile A=" <<AP
267 <<", Z=" <<ZP
268 <<", radius = " <<rP/fermi <<" fm"
269 <<G4endl;
270 G4cout <<"Initial target A=" <<AT
271 <<", Z=" <<ZT
272 <<", radius = " <<rT/fermi <<" fm"
273 <<G4endl;
274 G4cout <<"Projectile momentum and Energy/n
275 }
276 //
277 //
278 // The following variables are used to determi
279 // near-field (i.e. taking into consideration
280 //
281 G4double rm = ZP * ZT * elm_coupling / (E
282 G4double r = 0.0;
283 G4double rsq = 0.0;
284 //
285 //
286 // Initialise some of the variables which wll
287 // length for nucleons in the projectile and t
288 // number of abraded nucleons and the excitati
289 //
290 G4NuclearAbrasionGeometry *theAbrasionGeomet
291 G4double CT = 0.0;
292 G4double F = 0.0;
293 G4int Dabr = 0;
294 //
295 //
296 // The following loop is performed until the n
297 // abraded by the process is >1, i.e. an inter
298 //
299 G4bool skipInteraction = false; // It will
300 const G4int maxNumberOfLoops = 1000;
301 G4int loopCounter = -1;
302 while (Dabr == 0 && ++loopCounter < maxNumbe
303 {
304 //
305 //
306 // Sample the impact parameter. For the momen
307 // electrostatic effects on the impact paramet
308 // does not make any correction for the effect
309 //
310 G4double rPT = rP + rT;
311 G4double rPTsq = rPT * rPT;
312 //
313 //
314 // This is a "catch" to make sure we don't go
315 // energy is too low to overcome nuclear repul
316 // value of rm < fradius * rPT then we're unli
317 // impact parameter (energy of incident partic
318 //
319 if (rm >= fradius * rPT) {
320 skipInteraction = true;
321 }
322 //
323 //
324 // Now sample impact parameter until the crite
325 // and target overlap, but repulsion is taken
326 //
327 G4int evtcnt = 0;
328 r = 1.1 * rPT;
329 while (r > rPT && ++evtcnt < 1000) /* Loo
330 {
331 G4double bsq = rPTsq * G4UniformRand();
332 r = (rm + std::sqrt(rm*rm + 4
333 }
334 //
335 //
336 // We've tried to sample this 1000 times, but
337 //
338 if (evtcnt >= 1000) {
339 skipInteraction = true;
340 }
341
342 rsq = r * r;
343 //
344 //
345 // Now determine the chord-length through the
346 //
347 if (rT > rP)
348 {
349 G4double x = (rPsq + rsq - rTsq) / 2.0 /
350 if (x > 0.0) CT = 2.0 * std::sqrt(rTsq -
351 else CT = 2.0 * std::sqrt(rTsq -
352 }
353 else
354 {
355 G4double x = (rTsq + rsq - rPsq) / 2.0 /
356 if (x > 0.0) CT = 2.0 * std::sqrt(rTsq -
357 else CT = 2.0 * rT;
358 }
359 //
360 //
361 // Determine the number of abraded nucleons.
362 // abraded nucleons is used to sample the Pois
363 // distribution is sampled only ten times with
364 // and if it fails after this to find a case f
365 // nucleons >1, the impact parameter is re-sam
366 //
367 delete theAbrasionGeometry;
368 theAbrasionGeometry = new G4NuclearAbrasio
369 F = theAbrasionGeometry-
370 G4double lambda = 16.6*fermi / G4Pow::
371 G4double Mabr = F * AP * (1.0 - G4Ex
372 G4long n = 0;
373 for (G4int i = 0; i<10; ++i)
374 {
375 n = G4Poisson(Mabr);
376 if (n > 0)
377 {
378 if (n>AP) Dabr = (G4int) AP;
379 else Dabr = (G4int) n;
380 break;
381 }
382 }
383 } // End of while loop
384
385 if ( loopCounter >= maxNumberOfLoops || skip
386 // Assume nuclei do not collide and return
387 theParticleChange.SetStatusChange(isAlive)
388 theParticleChange.SetEnergyChange(theTrack
389 theParticleChange.SetMomentumChange(theTra
390 if (verboseLevel >= 2) {
391 G4cout <<"Particle energy too low to ove
392 G4cout <<"Event rejected and original tr
393 G4cout <<"##############################
394 <<"##############################
395 <<G4endl;
396 }
397 delete theAbrasionGeometry;
398 return &theParticleChange;
399 }
400
401 if (verboseLevel >= 2)
402 {
403 G4cout <<G4endl;
404 G4cout <<"Impact parameter = " <<r/ferm
405 G4cout <<"# Abraded nucleons = " <<Dabr <
406 }
407 //
408 //
409 // The number of abraded nucleons must be no g
410 // nucleons in either the projectile or the ta
411 // AT - Dabr < 2 then either we have only a nu
412 // projectile/target or we've tried to abrade
413 // should be limited.
414 //
415 if (AP - (G4double) Dabr < 2.0) Dabr = (G4in
416 if (AT - (G4double) Dabr < 2.0) Dabr = (G4in
417 //
418 //
419 // Determine the abraded secondary nucleons fr
420 // is a pointer to the prefragment from the pr
421 // of nucleons in theParticleChange which have
422 // from these is determined.
423 //
424 G4ThreeVector boost = pP.findBoostToCM();
425 G4Fragment *fragmentP = GetAbradedNucleons (
426 G4int nSecP = (G4int)theParticleCh
427 G4int i = 0;
428 for (i=0; i<nSecP; ++i)
429 {
430 TotalEPost += theParticleChange.GetSeconda
431 GetParticle()->GetTotalEnergy();
432 }
433 //
434 //
435 // Determine the number of spectators in the i
436 // projectile.
437 //
438 G4int DspcP = (G4int) (AP*F) - Dabr;
439 if (DspcP <= 0) DspcP = 0;
440 else if (DspcP > AP-Dabr) DspcP = ((G4int) A
441 //
442 //
443 // Determine excitation energy associated with
444 // projectile (EsP) and the excitation due to
445 // retained within the projectile (ExP). Add
446 // nucleus to the total energy of the secondar
447 //
448 G4bool excitationAbsorbedByProjectile = fals
449 if (fragmentP != nullptr)
450 {
451 G4double EsP = theAbrasionGeometry->GetExc
452 G4double ExP = 0.0;
453 if (Dabr < AT)
454 excitationAbsorbedByProjectile = G4Unifo
455 if (excitationAbsorbedByProjectile)
456 ExP = GetNucleonInducedExcitation(rP, rT
457 G4double xP = EsP + ExP;
458 if (xP > B*(AP-Dabr)) xP = B*(AP-Dabr);
459 G4LorentzVector lorentzVector = fragmentP-
460 lorentzVector.setE(lorentzVector.e()+xP);
461 fragmentP->SetMomentum(lorentzVector);
462 TotalEPost += lorentzVector.e();
463 }
464 G4double EMassP = TotalEPost;
465 //
466 //
467 // Determine the abraded secondary nucleons fr
468 // assumed that the same number of nucleons ar
469 // the projectile, and obviously no boost is a
470 // is a pointer to the prefragment from the ta
471 // of nucleons in theParticleChange which have
472 // from these is determined.
473 //
474 G4Fragment *fragmentT = GetAbradedNucleons (
475 G4int nSec = (G4int)theParticleChange.GetNum
476 for (i=nSecP; i<nSec; ++i)
477 {
478 TotalEPost += theParticleChange.GetSeconda
479 GetParticle()->GetTotalEnergy();
480 }
481 //
482 //
483 // Determine the number of spectators in the i
484 // target.
485 //
486 G4int DspcT = (G4int) (AT*F) - Dabr;
487 if (DspcT <= 0) DspcT = 0;
488 else if (DspcT > AP-Dabr) DspcT = ((G4int) A
489 //
490 //
491 // Determine excitation energy associated with
492 // target (EsT) and the excitation due to scat
493 // retained within the target (ExT). Add the
494 // nucleus to the total energy of the secondar
495 //
496 if (fragmentT != nullptr)
497 {
498 G4double EsT = theAbrasionGeometry->GetExc
499 G4double ExT = 0.0;
500 if (!excitationAbsorbedByProjectile)
501 ExT = GetNucleonInducedExcitation(rT, rP
502 G4double xT = EsT + ExT;
503 if (xT > B*(AT-Dabr)) xT = B*(AT-Dabr);
504 G4LorentzVector lorentzVector = fragmentT-
505 lorentzVector.setE(lorentzVector.e()+xT);
506 fragmentT->SetMomentum(lorentzVector);
507 TotalEPost += lorentzVector.e();
508 }
509 //
510 //
511 // Now determine the difference between the pr
512 // energy - this will be used to determine the
513 // of energy is to be imposed/attempted.
514 //
515 G4double deltaE = TotalEPre - TotalEPost;
516 if (deltaE > 0.0 && conserveEnergy)
517 {
518 G4double beta = std::sqrt(1.0 - EMassP*EMa
519 boost = boost / boost.mag() * beta;
520 }
521 //
522 //
523 // Now boost the secondaries from the projecti
524 //
525 G4ThreeVector pBalance = pP.vect();
526 for (i=0; i<nSecP; ++i)
527 {
528 G4DynamicParticle *dynamicP = theParticleC
529 GetParticle();
530 G4LorentzVector lorentzVector = dynamicP->
531 lorentzVector.boost(-boost);
532 dynamicP->Set4Momentum(lorentzVector);
533 pBalance -= lorentzVector.vect();
534 }
535 //
536 //
537 // Set the boost for the projectile prefragmen
538 // conservation of momentum. However, if the
539 // prefragment is not to be conserved this sim
540 // original projectile times the ratio of the
541 // of the prefragment (the excitation increase
542 // prefragment, and therefore modifying the bo
543 // the momentum of the prefragment being exces
544 //
545 if (fragmentP != nullptr)
546 {
547 G4LorentzVector lorentzVector = fragmentP-
548 G4double fragmentM = lorentzVec
549 if (conserveMomentum)
550 fragmentP->SetMomentum
551 (G4LorentzVector(pBalance,std::sqrt(pB
552 else
553 {
554 G4double fragmentGroundStateM = fragment
555 fragmentP->SetMomentum(lorentzVector.boo
556 }
557 }
558 //
559 //
560 // Output information to user if verbose infor
561 //
562 if (verboseLevel >= 2)
563 {
564 G4cout <<G4endl;
565 G4cout <<"--------------------------------
566 G4cout <<"Secondary nucleons from projecti
567 G4cout <<"--------------------------------
568 G4cout.precision(7);
569 for (i=0; i<nSecP; ++i)
570 {
571 G4cout <<"Particle # " <<i <<G4endl;
572 theParticleChange.GetSecondary(i)->GetPa
573 G4DynamicParticle *dyn = theParticleChan
574 G4cout <<"New nucleon (P) " <<dyn->GetDe
575 <<" : " <<dyn->Get4M
576 <<G4endl;
577 }
578 G4cout <<"---------------------------" <<G
579 G4cout <<"The projectile prefragment:" <<G
580 G4cout <<"---------------------------" <<G
581 if (fragmentP != nullptr)
582 G4cout <<*fragmentP <<G4endl;
583 else
584 G4cout <<"(No residual prefragment)" <<G
585 G4cout <<G4endl;
586 G4cout <<"-------------------------------"
587 G4cout <<"Secondary nucleons from target:"
588 G4cout <<"-------------------------------"
589 G4cout.precision(7);
590 for (i=nSecP; i<nSec; ++i)
591 {
592 G4cout <<"Particle # " <<i <<G4endl;
593 theParticleChange.GetSecondary(i)->GetPa
594 G4DynamicParticle *dyn = theParticleChan
595 G4cout <<"New nucleon (T) " <<dyn->GetDe
596 <<" : " <<dyn->Get4M
597 <<G4endl;
598 }
599 G4cout <<"-----------------------" <<G4end
600 G4cout <<"The target prefragment:" <<G4end
601 G4cout <<"-----------------------" <<G4end
602 if (fragmentT != nullptr)
603 G4cout <<*fragmentT <<G4endl;
604 else
605 G4cout <<"(No residual prefragment)" <<G
606 }
607 //
608 //
609 // Now we can decay the nuclear fragments if p
610 // collected and boosted as well. This is per
611 //
612 if (fragmentP !=nullptr)
613 {
614 G4ReactionProductVector *products = nullpt
615 // if (fragmentP->GetZ_asInt() != fragm
616 products = theExcitationHandler->BreakItUp
617 // else
618 // products = theExcitationHandlerx->Bre
619 delete fragmentP;
620 fragmentP = nullptr;
621
622 G4ReactionProductVector::iterator iter;
623 for (iter = products->begin(); iter != pro
624 {
625 G4DynamicParticle *secondary =
626 new G4DynamicParticle((*iter)->GetDefi
627 (*iter)->GetTotalEnergy(), (*iter)->Ge
628 theParticleChange.AddSecondary (secondar
629 G4String particleName = (*iter)->GetDefi
630 delete (*iter); // get rid of leftover p
631 if (verboseLevel >= 2 && particleName.fi
632 {
633 G4cout <<"------------------------" <<
634 G4cout <<"The projectile fragment:" <<
635 G4cout <<"------------------------" <<
636 G4cout <<" fragmentP = " <<particleNam
637 <<" Energy = " <<secondary->
638 <<G4endl;
639 }
640 }
641 delete products;
642 }
643 //
644 //
645 // Now decay the target nucleus - no boost is
646 // approximation it is assumed that there is n
647 // the projectile.
648 //
649 if (fragmentT != nullptr)
650 {
651 G4ReactionProductVector *products = nullpt
652 // if (fragmentT->GetZ_asInt() != fragm
653 products = theExcitationHandler->BreakIt
654 // else
655 // products = theExcitationHandlerx->Bre
656 // delete fragmentT;
657 fragmentT = nullptr;
658
659 G4ReactionProductVector::iterator iter;
660 for (iter = products->begin(); iter != pro
661 {
662 G4DynamicParticle *secondary =
663 new G4DynamicParticle((*iter)->GetDefi
664 (*iter)->GetTotalEnergy(), (*iter)->Ge
665 theParticleChange.AddSecondary (secondar
666 G4String particleName = (*iter)->GetDefi
667 delete (*iter); // get rid of leftover p
668 if (verboseLevel >= 2 && particleName.fi
669 {
670 G4cout <<"--------------------" <<G4en
671 G4cout <<"The target fragment:" <<G4en
672 G4cout <<"--------------------" <<G4en
673 G4cout <<" fragmentT = " <<particleNam
674 <<" Energy = " <<secondary->
675 <<G4endl;
676 }
677 }
678 delete products;
679 }
680
681 if (verboseLevel >= 2)
682 G4cout <<"###############################
683 <<"###############################
684 <<G4endl;
685
686 delete theAbrasionGeometry;
687 return &theParticleChange;
688 }
689 //////////////////////////////////////////////
690 //
691 G4Fragment *G4WilsonAbrasionModel::GetAbradedN
692 G4double Z, G4double r)
693 {
694 //
695 //
696 // Initialise variables. tau is the Fermi rad
697 // p..., C... and gamma are used to help sampl
698 // spectrum.
699 //
700
701 G4double pK = hbarc * G4Pow::GetInstance()
702 if (A <= 24.0) pK *= -0.229*G4Pow::GetInstan
703 G4double pKsq = pK * pK;
704 G4double p1sq = 2.0/5.0 * pKsq;
705 G4double p2sq = 6.0/5.0 * pKsq;
706 G4double p3sq = 500.0 * 500.0;
707 G4double C1 = 1.0;
708 G4double C2 = 0.03;
709 G4double C3 = 0.0002;
710 G4double gamma = 90.0 * MeV;
711 G4double maxn = C1 + C2 + C3;
712 //
713 //
714 // initialise the number of secondary nucleons
715 // the type of nucleon abraded to proton ... j
716 //
717 G4double Aabr = 0.0;
718 G4double Zabr = 0.0;
719 G4ParticleDefinition *typeNucleon = G4Proton
720 G4DynamicParticle *dynamicNucleon = nullptr;
721 G4ParticleMomentum pabr(0.0, 0.0, 0.0);
722 //
723 //
724 // Now go through each abraded nucleon and sam
725 //
726 G4bool isForLoopExitAnticipated = false;
727 for (G4int i=0; i<Dabr; ++i)
728 {
729 //
730 //
731 // Sample the nucleon momentum distribution by
732 // reject values of p == 0.0 since this causes
733 //
734 G4double p = 0.0;
735 G4bool found = false;
736 const G4int maxNumberOfLoops = 100000;
737 G4int loopCounter = -1;
738 while (!found && ++loopCounter < maxNumber
739 {
740 while (p <= 0.0) p = npK * pK * G4Unifor
741 G4double psq = p * p;
742 found = maxn * G4UniformRand() < C1*G4Ex
743 C2*G4Exp(-psq/p2sq/2.0) + C3*G4Exp(-ps
744 }
745 if ( loopCounter >= maxNumberOfLoops )
746 {
747 isForLoopExitAnticipated = true;
748 break;
749 }
750 //
751 //
752 // Determine the type of particle abraded. Ca
753 // and the probability is determine to be prop
754 // in the nucleus at each stage.
755 //
756 G4double prob = (Z-Zabr)/(A-Aabr);
757 if (G4UniformRand()<prob)
758 {
759 Zabr++;
760 typeNucleon = G4Proton::ProtonDefinition
761 }
762 else
763 typeNucleon = G4Neutron::NeutronDefiniti
764 Aabr++;
765 //
766 //
767 // The angular distribution of the secondary n
768 // isotropic distribution in the rest frame of
769 // boosted later.
770 //
771 G4double costheta = 2.*G4UniformRand()-1.0
772 G4double sintheta = std::sqrt((1.0 - costh
773 G4double phi = 2.0*pi*G4UniformRand()
774 G4ThreeVector direction(sintheta*std::cos(
775 G4double nucleonMass = typeNucleon->GetPDG
776 G4double E = std::sqrt(p*p + nuc
777 dynamicNucleon = new G4DynamicParticle(typ
778 theParticleChange.AddSecondary (dynamicNuc
779 pabr += p*direction;
780 }
781 //
782 //
783 // Next determine the details of the nuclear p
784 // is one or more protons in the residue. (No
785 // energy is a safety factor to avoid any poss
786 // energy.)
787 //
788 G4Fragment *fragment = nullptr;
789 if ( ! isForLoopExitAnticipated && Z-Zabr>=1
790 {
791 G4double ionMass = G4ParticleTable::GetPar
792 GetIonMass(G4lrint(Z-Za
793 G4double E = std::sqrt(pabr.mag2() +
794 G4LorentzVector lorentzVector = G4LorentzV
795 fragment =
796 new G4Fragment((G4int) (A-Aabr), (G4int)
797 }
798
799 return fragment;
800 }
801 //////////////////////////////////////////////
802 //
803 G4double G4WilsonAbrasionModel::GetNucleonIndu
804 (G4double rP, G4double rT, G4double r)
805 {
806 //
807 //
808 // Initialise variables.
809 //
810 G4double Cl = 0.0;
811 G4double rPsq = rP * rP;
812 G4double rTsq = rT * rT;
813 G4double rsq = r * r;
814 //
815 //
816 // Depending upon the impact parameter, a diff
817 // is used.
818 //
819 if (r > rT) Cl = 2.0*std::sqrt(rPsq + 2.0*r*
820 else Cl = 2.0*rP;
821 //
822 //
823 // The next lines have been changed to include
824 // projectile and target are too close, Ct is
825 // Otherwise the standard Wilson algorithm sho
826 // PRT 20091023.
827 //
828 G4double Ct = 0.0;
829 if (rT > rP && rsq < rTsq - rPsq) Ct =
830 else if (rP > rT && rsq < rPsq - rTsq) Ct =
831 else {
832 G4double bP = (rPsq+rsq-rTsq)/2.0/r;
833 G4double x = rPsq - bP*bP;
834 if (x < 0.0) {
835 G4cerr <<"##############################
836 <<"##############################
837 <<G4endl;
838 G4cerr <<"ERROR IN G4WilsonAbrasionModel
839 <<G4endl;
840 G4cerr <<"rPsq - bP*bP < 0.0 and cannot
841 G4cerr <<"Set to zero instead" <<G4endl;
842 G4cerr <<"##############################
843 <<"##############################
844 <<G4endl;
845 }
846 Ct = 2.0*std::sqrt(x);
847 }
848
849 G4double Ex = 13.0 * Cl / fermi;
850 if (Ct > 1.5*fermi)
851 Ex += 13.0 * Cl / fermi /3.0 * (Ct/fermi -
852
853 return Ex;
854 }
855 //////////////////////////////////////////////
856 //
857 void G4WilsonAbrasionModel::SetUseAblation (G4
858 {
859 if (useAblation != useAblation1)
860 {
861 useAblation = useAblation1;
862 if (useAblation)
863 {
864 theAblation = new G4WilsonAblationModel;
865 theAblation->SetVerboseLevel(verboseLeve
866 theExcitationHandler->SetEvaporation(the
867 }
868 else
869 {
870 delete theExcitationHandler;
871 theAblation = nullp
872 theExcitationHandler = new G4Excitation
873 }
874 }
875 return;
876 }
877 //////////////////////////////////////////////
878 //
879 void G4WilsonAbrasionModel::PrintWelcomeMessag
880 {
881 G4cout <<G4endl;
882 G4cout <<" *********************************
883 <<G4endl;
884 G4cout <<" Nuclear abrasion model for nuclea
885 <<G4endl;
886 G4cout <<" (Written by QinetiQ Ltd for the E
887 <<G4endl;
888 G4cout <<" *********************************
889 <<G4endl;
890 G4cout << G4endl;
891
892 return;
893 }
894 //////////////////////////////////////////////
895 //
896