| 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: G4WilsonAblationModel.
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 // 08 December 2009, P R Truscott, QinetiQ Ltd
59 // Ver 1.0
60 // Updated as a result of changes in the G4Eva
61 // affect mostly SelectSecondariesByEvaporatio
62 // associated with the evaporation model which
63 // OPTxs to select the inverse cross-sectio
64 // OPTxs = 0 => Dostrovski's parameter
65 // OPTxs = 1 or 2 => Chatterjee's paramater
66 // OPTxs = 3 or 4 => Kalbach's parameteriza
67 // useSICB => use superimposed Coulo
68 // sections
69 // Other problem found with G4Fragment definit
70 // **G4ParticleDefinition**. This does not al
71 // fragment for some reason. Now the fragment
72 // G4Fragment *fragment = new G4Fragment(A,
73 // to avoid this quirk. Bug found in SelectSe
74 // equated to evapType[i] whereas previously i
75 //
76 // 06 August 2015, A. Ribon, CERN
77 // Migrated std::exp and std::pow to the faste
78 //
79 // 09 June 2017, C. Mancini Terracciano, INFN
80 // Fixed bug on the initialization of Photon E
81 //
82 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
83 //////////////////////////////////////////////
84 //
85 #include <iomanip>
86 #include <numeric>
87
88 #include "G4WilsonAblationModel.hh"
89 #include "G4PhysicalConstants.hh"
90 #include "G4SystemOfUnits.hh"
91 #include "Randomize.hh"
92 #include "G4ParticleTable.hh"
93 #include "G4IonTable.hh"
94 #include "G4Alpha.hh"
95 #include "G4He3.hh"
96 #include "G4Triton.hh"
97 #include "G4Deuteron.hh"
98 #include "G4Proton.hh"
99 #include "G4Neutron.hh"
100 #include "G4AlphaEvaporationChannel.hh"
101 #include "G4He3EvaporationChannel.hh"
102 #include "G4TritonEvaporationChannel.hh"
103 #include "G4DeuteronEvaporationChannel.hh"
104 #include "G4ProtonEvaporationChannel.hh"
105 #include "G4NeutronEvaporationChannel.hh"
106 #include "G4PhotonEvaporation.hh"
107 #include "G4LorentzVector.hh"
108 #include "G4VEvaporationChannel.hh"
109
110 #include "G4Exp.hh"
111 #include "G4Pow.hh"
112
113 #include "G4PhysicsModelCatalog.hh"
114
115 //////////////////////////////////////////////
116 //
117 G4WilsonAblationModel::G4WilsonAblationModel()
118 {
119 //
120 //
121 // Send message to stdout to advise that the G
122 //
123 PrintWelcomeMessage();
124 //
125 //
126 // Set the default verbose level to 0 - no out
127 //
128 verboseLevel = 0;
129 //
130 //
131 // Set the binding energy per nucleon .... did
132 // model for nuclear de-excitation?
133 //
134 B = 10.0 * MeV;
135 //
136 //
137 // It is possuble to switch off secondary part
138 // final nuclear fragment). The default is on
139 //
140 produceSecondaries = true;
141 //
142 //
143 // Now we need to define the decay modes. We'
144 // to help determine the kinematics of the dec
145 //
146 nFragTypes = 6;
147 fragType[0] = G4Alpha::Alpha();
148 fragType[1] = G4He3::He3();
149 fragType[2] = G4Triton::Triton();
150 fragType[3] = G4Deuteron::Deuteron();
151 fragType[4] = G4Proton::Proton();
152 fragType[5] = G4Neutron::Neutron();
153 for(G4int i=0; i<200; ++i) { fSig[i] = 0.0;
154 //
155 //
156 // Set verboseLevel default to no output.
157 //
158 verboseLevel = 0;
159 theChannelFactory = new G4EvaporationFactory
160 theChannels = theChannelFactory->GetChannel(
161 //
162 //
163 // Set defaults for evaporation classes. Thes
164 // "set" methods.
165 //
166 OPTxs = 3;
167 useSICB = false;
168 fragmentVector = 0;
169
170 secID = G4PhysicsModelCatalog::GetModelID("m
171 }
172 //////////////////////////////////////////////
173 //
174 G4WilsonAblationModel::~G4WilsonAblationModel(
175 {}
176
177 //////////////////////////////////////////////
178 //
179 G4FragmentVector *G4WilsonAblationModel::Break
180 (const G4Fragment &theNucleus)
181 {
182 //
183 //
184 // Initilise the pointer to the G4FragmentVect
185 // about the breakup.
186 //
187 fragmentVector = new G4FragmentVector;
188 fragmentVector->clear();
189 //
190 //
191 // Get the A, Z and excitation of the nucleus.
192 //
193 G4int A = theNucleus.GetA_asInt();
194 G4int Z = theNucleus.GetZ_asInt();
195 G4double ex = theNucleus.GetExcitationEnergy
196 if (verboseLevel >= 2)
197 {
198 G4cout <<"oooooooooooooooooooooooooooooooo
199 <<"oooooooooooooooooooooooooooooooo
200 <<G4endl;
201 G4cout.precision(6);
202 G4cout <<"IN G4WilsonAblationModel" <<G4en
203 G4cout <<"Initial prefragment A=" <<A
204 <<", Z=" <<Z
205 <<", excitation energy = " <<ex/MeV
206 <<G4endl;
207 }
208 //
209 //
210 // Check that there is a nucleus to speak of.
211 // or its just a proton or neutron. In either
212 // (from the Lorentz vector) is not used.
213 //
214 if (A == 0)
215 {
216 if (verboseLevel >= 2)
217 {
218 G4cout <<"No nucleus to decay" <<G4endl;
219 G4cout <<"oooooooooooooooooooooooooooooo
220 <<"oooooooooooooooooooooooooooooo
221 <<G4endl;
222 }
223 return fragmentVector;
224 }
225 else if (A == 1)
226 {
227 G4LorentzVector lorentzVector = theNucleus
228 lorentzVector.setE(lorentzVector.e()-ex+10
229 if (Z == 0)
230 {
231 G4Fragment *fragment = new G4Fragment(lo
232 if (fragment != nullptr) { fragment->Set
233 fragmentVector->push_back(fragment);
234 }
235 else
236 {
237 G4Fragment *fragment = new G4Fragment(lo
238 if (fragment != nullptr) { fragment->Set
239 fragmentVector->push_back(fragment);
240 }
241 if (verboseLevel >= 2)
242 {
243 G4cout <<"Final fragment is in fact only
244 G4cout <<(*fragmentVector)[0] <<G4endl;
245 G4cout <<"oooooooooooooooooooooooooooooo
246 <<"oooooooooooooooooooooooooooooo
247 <<G4endl;
248 }
249 return fragmentVector;
250 }
251 //
252 //
253 // Then the number of nucleons ablated (either
254 // fragments) is based on a simple argument fo
255 //
256 G4int DAabl = (G4int) (ex / B);
257 if (DAabl > A) DAabl = A;
258 // The following lines are no longer accurate
259 // if (verboseLevel >= 2)
260 // G4cout <<"Number of nucleons ejected = "
261
262 //
263 //
264 // Determine the nuclear fragment from the abl
265 // Rudstam equation.
266 //
267 G4int AF = A - DAabl;
268 G4int ZF = 0;
269
270 if (AF > 0)
271 {
272 G4Pow* g4calc = G4Pow::GetInstance();
273 G4double AFd = (G4double) AF;
274 G4double R = 11.8 / g4calc->powZ(AF, 0.45)
275 G4int minZ = std::max(1, Z - DAabl);
276 //
277 //
278 // Here we define an integral probability dist
279 // equation assuming a constant AF.
280 //
281 G4int zmax = std::min(199, Z);
282 G4double sum = 0.0;
283 for (ZF=minZ; ZF<=zmax; ++ZF)
284 {
285 sum += G4Exp(-R*g4calc->powA(std::abs(ZF
286 fSig[ZF] = sum;
287 }
288 //
289 //
290 // Now sample that distribution to determine a
291 //
292 sum *= G4UniformRand();
293 for (ZF=minZ; ZF<=zmax; ++ZF) {
294 if(sum <= fSig[ZF]) { break; }
295 }
296 }
297 G4int DZabl = Z - ZF;
298 //
299 //
300 // Now determine the nucleons or nuclei which
301 // is for the production of alphas, then other
302 // binding energy. The energies assigned to th
303 // provisional for the moment (the 10eV is jus
304 // excitation energies due to rounding).
305 //
306 G4double totalEpost = 0.0;
307 evapType.clear();
308 for (G4int ift=0; ift<nFragTypes; ift++)
309 {
310 G4ParticleDefinition *type = fragType[ift]
311 G4double n = std::floor((G4double) DAabl
312 G4double n1 = 1.0E+10;
313 if (fragType[ift]->GetPDGCharge() > 0.0)
314 n1 = std::floor((G4double) DZabl / type-
315 if (n > n1) n = n1;
316 if (n > 0.0)
317 {
318 G4double mass = type->GetPDGMass();
319 for (G4int j=0; j<(G4int) n; j++)
320 {
321 totalEpost += mass;
322 evapType.push_back(type);
323 }
324 DAabl -= (G4int) (n * type->GetBaryonNum
325 DZabl -= (G4int) (n * type->GetPDGCharge
326 }
327 }
328 //
329 //
330 // Determine the properties of the final nucle
331 // the final fragment is predicted to have a n
332 // really it's the particle last in the vector
333 // final fragment. Therefore delete this from
334 // case.
335 //
336 G4double massFinalFrag = 0.0;
337 if (AF > 0)
338 massFinalFrag = G4ParticleTable::GetPartic
339 GetIonMass(ZF,AF);
340 else
341 {
342 G4ParticleDefinition *type = evapType[evap
343 AF = type->GetBary
344 ZF = (G4int) (type
345 evapType.erase(evapType.end()-1);
346 }
347 totalEpost += massFinalFrag;
348 //
349 //
350 // Provide verbose output on the nuclear fragm
351 //
352 if (verboseLevel >= 2)
353 {
354 G4cout <<"Final fragment A=" <<AF
355 <<", Z=" <<ZF
356 <<G4endl;
357 for (G4int ift=0; ift<nFragTypes; ift++)
358 {
359 G4ParticleDefinition *type = fragType[if
360 G4long n = std::count(evapType.cbegin(),
361 if (n > 0)
362 G4cout <<"Particle type: " <<std::setw
363 <<", number of particles emitte
364 }
365 }
366 //
367 // Add the total energy from the fragment. No
368 // to be de-excited and does not undergo photo
369 // this is a bit of a crude model?
370 //
371 G4double massPreFrag = theNucleus.GetGr
372 G4double totalEpre = massPreFrag + ex
373 G4double excess = totalEpre - tota
374 // G4Fragment *resultNucleus(theNucleus);
375 G4Fragment *resultNucleus = new G4Fragment(A
376 G4ThreeVector boost(0.0,0.0,0.0);
377 std::size_t nEvap = 0;
378 if (produceSecondaries && evapType.size()>0)
379 {
380 if (excess > 0.0)
381 {
382 SelectSecondariesByEvaporation (resultNu
383 nEvap = fragmentVector->size();
384 boost = resultNucleus->GetMomentum().fin
385 if (evapType.size() > 0)
386 SelectSecondariesByDefault (boost);
387 }
388 else
389 SelectSecondariesByDefault(G4ThreeVector
390 }
391
392 if (AF > 0)
393 {
394 G4double mass = G4ParticleTable::GetPartic
395 GetIonMass(ZF,AF);
396 G4double e = mass + 10.0*eV;
397 G4double p = std::sqrt(e*e-mass*mass);
398 G4ThreeVector direction(0.0,0.0,1.0);
399 G4LorentzVector lorentzVector = G4LorentzV
400 lorentzVector.boost(-boost);
401 G4Fragment* frag = new G4Fragment(AF, ZF,
402 if (frag != nullptr) { frag->SetCreatorMod
403 fragmentVector->push_back(frag);
404 }
405 delete resultNucleus;
406 //
407 //
408 // Provide verbose output on the ablation prod
409 //
410 if (verboseLevel >= 2)
411 {
412 if (nEvap > 0)
413 {
414 G4cout <<"----------------------" <<G4en
415 G4cout <<"Evaporated particles :" <<G4en
416 G4cout <<"----------------------" <<G4en
417 }
418 std::size_t ie = 0;
419 for (auto iter = fragmentVector->cbegin();
420 iter != fragmentVector->cend();
421 {
422 if (ie == nEvap)
423 {
424 // G4cout <<*iter <<G4endl;
425 G4cout <<"----------------------------
426 G4cout <<"Particles from default emiss
427 G4cout <<"----------------------------
428 }
429 G4cout <<*iter <<G4endl;
430 }
431 G4cout <<"oooooooooooooooooooooooooooooooo
432 <<"oooooooooooooooooooooooooooooooo
433 <<G4endl;
434 }
435
436 return fragmentVector;
437 }
438 //////////////////////////////////////////////
439 //
440 void G4WilsonAblationModel::SelectSecondariesB
441 (G4Fragment *intermediateNucleus)
442 {
443 G4Fragment theResidualNucleus = *intermediat
444 G4bool evaporate = true;
445 // Loop checking, 05-Aug-2015, Vladimir Ivan
446 while (evaporate && evapType.size() != 0)
447 {
448 //
449 //
450 // Here's the cheaky bit. We're hijacking the
451 // more accurately sample to kinematics, but t
452 // fragments will be the ones of our choosing
453 //
454 std::vector <G4VEvaporationChannel*> theC
455 theChannels1.clear();
456 std::vector <G4VEvaporationChannel*>::iter
457 VectorOfFragmentTypes::iterator iter;
458 std::vector <VectorOfFragmentTypes::iterat
459 iters.clear();
460 iter = std::find(evapType.begin(), evapTyp
461 if (iter != evapType.end())
462 {
463 theChannels1.push_back(new G4AlphaEvapor
464 i = theChannels1.end() - 1;
465 (*i)->SetOPTxs(OPTxs);
466 (*i)->UseSICB(useSICB);
467 // (*i)->Initialize(theResidualNucleus);
468 iters.push_back(iter);
469 }
470 iter = std::find(evapType.begin(), evapTyp
471 if (iter != evapType.end())
472 {
473 theChannels1.push_back(new G4He3Evaporat
474 i = theChannels1.end() - 1;
475 (*i)->SetOPTxs(OPTxs);
476 (*i)->UseSICB(useSICB);
477 // (*i)->Initialize(theResidualNucleus);
478 iters.push_back(iter);
479 }
480 iter = std::find(evapType.begin(), evapTyp
481 if (iter != evapType.end())
482 {
483 theChannels1.push_back(new G4TritonEvapo
484 i = theChannels1.end() - 1;
485 (*i)->SetOPTxs(OPTxs);
486 (*i)->UseSICB(useSICB);
487 // (*i)->Initialize(theResidualNucleus);
488 iters.push_back(iter);
489 }
490 iter = std::find(evapType.begin(), evapTyp
491 if (iter != evapType.end())
492 {
493 theChannels1.push_back(new G4DeuteronEva
494 i = theChannels1.end() - 1;
495 (*i)->SetOPTxs(OPTxs);
496 (*i)->UseSICB(useSICB);
497 // (*i)->Initialize(theResidualNucleus);
498 iters.push_back(iter);
499 }
500 iter = std::find(evapType.begin(), evapTyp
501 if (iter != evapType.end())
502 {
503 theChannels1.push_back(new G4ProtonEvapo
504 i = theChannels1.end() - 1;
505 (*i)->SetOPTxs(OPTxs);
506 (*i)->UseSICB(useSICB);
507 // (*i)->Initialize(theResidualNucleus);
508 iters.push_back(iter);
509 }
510 iter = std::find(evapType.begin(), evapTyp
511 if (iter != evapType.end())
512 {
513 theChannels1.push_back(new G4NeutronEvap
514 i = theChannels1.end() - 1;
515 (*i)->SetOPTxs(OPTxs);
516 (*i)->UseSICB(useSICB);
517 // (*i)->Initialize(theResidualNucleus);
518 iters.push_back(iter);
519 }
520 std::size_t nChannels = theChannels1.size(
521
522 G4double totalProb = 0.0;
523 G4int ich = 0;
524 G4double probEvapType[6] = {0.0};
525 for (auto iterEv=theChannels1.cbegin();
526 iterEv!=theChannels1.cend(); ++i
527 totalProb += (*iterEv)->GetEmissionProba
528 probEvapType[ich] = totalProb;
529 ++ich;
530 }
531 if (totalProb > 0.0) {
532 //
533 //
534 // The emission probability for at least one o
535 // positive, therefore work out which one shou
536 // the nucleus.
537 //
538 G4double xi = totalProb*G4UniformRand();
539 std::size_t ii = 0;
540 for (ii=0; ii<nChannels; ++ii)
541 {
542 if (xi < probEvapType[ii]) { break; }
543 }
544 if (ii >= nChannels) { ii = nChannels -
545 G4FragmentVector *evaporationResult = th
546 BreakUpFragment(intermediateNucleus);
547 if ((*evaporationResult)[0] != nullptr)
548 {
549 (*evaporationResult)[0]->SetCreatorMod
550 }
551 fragmentVector->push_back((*evaporationR
552 intermediateNucleus = (*evaporationResul
553 delete evaporationResult;
554 }
555 else
556 {
557 //
558 //
559 // Probability for further evaporation is nil
560 // routine and set the energies of the seconda
561 //
562 evaporate = false;
563 }
564 }
565
566 return;
567 }
568 //////////////////////////////////////////////
569 //
570 void G4WilsonAblationModel::SelectSecondariesB
571 {
572 for (std::size_t i=0; i<evapType.size(); ++i
573 {
574 G4ParticleDefinition *type = evapType[i];
575 G4double mass = type->GetPDGM
576 G4double e = mass + 10.0*e
577 G4double p = std::sqrt(e*e
578 G4double costheta = 2.0*G4Uniform
579 G4double sintheta = std::sqrt((1.
580 G4double phi = twopi * G4Uni
581 G4ThreeVector direction(sintheta*std::cos(
582 G4LorentzVector lorentzVector = G4LorentzV
583 lorentzVector.boost(-boost);
584 // Possibility that the following line is not
585 // from particle definition. Force values. P
586 // G4Fragment *fragment =
587 // new G4Fragment(lorentzVector, type);
588 G4int A = type->GetBaryonNumber();
589 G4int Z = (G4int) (type->GetPDGCharge() +
590 G4Fragment *fragment =
591 new G4Fragment(A, Z, lorentzVector);
592 if (fragment != nullptr) { fragment->SetCr
593 fragmentVector->push_back(fragment);
594 }
595 }
596 //////////////////////////////////////////////
597 //
598 void G4WilsonAblationModel::PrintWelcomeMessag
599 {
600 G4cout <<G4endl;
601 G4cout <<" *********************************
602 <<G4endl;
603 G4cout <<" Nuclear ablation model for nuclea
604 <<G4endl;
605 G4cout <<" (Written by QinetiQ Ltd for the E
606 <<G4endl;
607 G4cout <<" !!! WARNING: This model is not we
608 <<G4endl;
609 G4cout <<" *********************************
610 <<G4endl;
611 G4cout << G4endl;
612
613 return;
614 }
615 //////////////////////////////////////////////
616 //
617