| 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 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // Hadronic Process: Nuclear De-excitations
27 // by V. Lara (May 1998)
28 //
29 //
30 // Modified:
31 // 30 June 1998 by V. Lara:
32 // -Modified the Transform method for use
33 // therefore G4IonTable. It makes possib
34 // of fragments (G4Fragment) produced in
35 // G4DynamicParticle
36 // -It uses default algorithms for:
37 // Evaporation: G4Evaporation
38 // MultiFragmentation: G4StatMF
39 // Fermi Breakup model: G4FermiBr
40 // 24 Jul 2008 by M. A. Cortes Giraldo:
41 // -Max Z,A for Fermi Break-Up turns to 9
42 // -BreakItUp() reorganised and bug in Ev
43 // -Transform() optimised
44 // (September 2008) by J. M. Quesada. External
45 // -inverse cross section option (default
46 // -superimposed Coulomb barrier (if useS
47 // September 2009 by J. M. Quesada:
48 // -according to Igor Pshenichnov, SMM wi
49 // 27 Nov 2009 by V.Ivanchenko:
50 // -cleanup the logic, reduce number inte
51 // 11 May 2010 by V.Ivanchenko:
52 // -FermiBreakUp activated, used integer
53 // final photon deexcitation; used check
54 // unstable fragments with A <5
55 // 22 March 2011 by V.Ivanchenko: general clea
56 // products of Fermi Break Up cannot be
57 // 30 March 2011 by V.Ivanchenko removed priva
58 // to the source
59 // 23 January 2012 by V.Ivanchenko general cle
60 // objects, propagate G4PhotonEvaporation p
61 // not delete it here
62
63 #include "G4ExcitationHandler.hh"
64 #include "G4SystemOfUnits.hh"
65 #include "G4LorentzVector.hh"
66 #include "G4ThreeVector.hh"
67 #include "G4ParticleTable.hh"
68 #include "G4ParticleTypes.hh"
69 #include "G4Ions.hh"
70 #include "G4Electron.hh"
71 #include "G4Lambda.hh"
72
73 #include "G4VMultiFragmentation.hh"
74 #include "G4VFermiBreakUp.hh"
75 #include "G4Element.hh"
76 #include "G4ElementTable.hh"
77
78 #include "G4VEvaporation.hh"
79 #include "G4VEvaporationChannel.hh"
80 #include "G4Evaporation.hh"
81 #include "G4PhotonEvaporation.hh"
82 #include "G4StatMF.hh"
83 #include "G4FermiBreakUpVI.hh"
84 #include "G4NuclearLevelData.hh"
85 #include "G4PhysicsModelCatalog.hh"
86
87 G4ExcitationHandler::G4ExcitationHandler()
88 : minEForMultiFrag(1.*CLHEP::TeV), minExcita
89 maxExcitation(100.*CLHEP::MeV)
90 {
91 thePartTable = G4ParticleTable::GetParticleT
92 theTableOfIons = thePartTable->GetIonTable()
93 nist = G4NistManager::Instance();
94
95 theMultiFragmentation = new G4StatMF();
96 theFermiModel = new G4FermiBreakUpVI();
97 thePhotonEvaporation = new G4PhotonEvaporati
98 SetEvaporation(new G4Evaporation(thePhotonEv
99 theResults.reserve(60);
100 results.reserve(30);
101 theEvapList.reserve(30);
102
103 theElectron = G4Electron::Electron();
104 theNeutron = G4Neutron::NeutronDefinition();
105 theProton = G4Proton::ProtonDefinition();
106 theDeuteron = G4Deuteron::DeuteronDefinition
107 theTriton = G4Triton::TritonDefinition();
108 theHe3 = G4He3::He3Definition();
109 theAlpha = G4Alpha::AlphaDefinition();
110 theLambda = G4Lambda::Lambda();
111
112 fLambdaMass = theLambda->GetPDGMass();
113
114 if(fVerbose > 1) { G4cout << "### New handle
115 }
116
117 G4ExcitationHandler::~G4ExcitationHandler()
118 {
119 delete theMultiFragmentation;
120 delete theFermiModel;
121 if(isEvapLocal) { delete theEvaporation; }
122 }
123
124 void G4ExcitationHandler::SetParameters()
125 {
126 G4NuclearLevelData* ndata = G4NuclearLevelDa
127 auto param = ndata->GetParameters();
128 isActive = true;
129 // check if de-excitation is needed
130 if (fDummy == param->GetDeexChannelsType())
131 isActive = false;
132 } else {
133 // upload data for elements used in geomet
134 G4int Zmax = 20;
135 const G4ElementTable* table = G4Element::G
136 for (auto const & elm : *table) { Zmax = s
137 ndata->UploadNuclearLevelData(Zmax+1);
138 }
139 minEForMultiFrag = param->GetMinExPerNucleou
140 minExcitation = param->GetMinExcitation();
141 maxExcitation = param->GetPrecoHighEnergy();
142
143 // allowing local debug printout
144 fVerbose = std::max(fVerbose, param->GetVerb
145 if (isActive) {
146 if (nullptr == thePhotonEvaporation) {
147 SetPhotonEvaporation(new G4PhotonEvapora
148 }
149 if (nullptr == theFermiModel) {
150 SetFermiModel(new G4FermiBreakUpVI());
151 }
152 if (nullptr == theMultiFragmentation) {
153 SetMultiFragmentation(new G4StatMF());
154 }
155 if (nullptr == theEvaporation) {
156 SetEvaporation(new G4Evaporation(thePhot
157 }
158 }
159 theFermiModel->SetVerbose(fVerbose);
160 if(fVerbose > 1) {
161 G4cout << "G4ExcitationHandler::SetParamet
162 }
163 }
164
165 void G4ExcitationHandler::Initialise()
166 {
167 if(isInitialised) { return; }
168 if(fVerbose > 1) {
169 G4cout << "G4ExcitationHandler::Initialise
170 }
171 G4DeexPrecoParameters* param =
172 G4NuclearLevelData::GetInstance()->GetPara
173 isInitialised = true;
174 SetParameters();
175 if(isActive) {
176 theFermiModel->Initialise();
177 theEvaporation->InitialiseChannels();
178 }
179 // dump level is controlled by parameter cla
180 param->Dump();
181 }
182
183 void G4ExcitationHandler::SetEvaporation(G4VEv
184 {
185 if(nullptr != ptr && ptr != theEvaporation)
186 theEvaporation = ptr;
187 theEvaporation->SetPhotonEvaporation(thePh
188 theEvaporation->SetFermiBreakUp(theFermiMo
189 isEvapLocal = flag;
190 if(fVerbose > 1) {
191 G4cout << "G4ExcitationHandler::SetEvapo
192 }
193 }
194 }
195
196 void
197 G4ExcitationHandler::SetMultiFragmentation(G4V
198 {
199 if(nullptr != ptr && ptr != theMultiFragment
200 delete theMultiFragmentation;
201 theMultiFragmentation = ptr;
202 }
203 }
204
205 void G4ExcitationHandler::SetFermiModel(G4VFer
206 {
207 if(nullptr != ptr && ptr != theFermiModel) {
208 delete theFermiModel;
209 theFermiModel = ptr;
210 if(nullptr != theEvaporation) {
211 theEvaporation->SetFermiBreakUp(theFermi
212 }
213 }
214 }
215
216 void
217 G4ExcitationHandler::SetPhotonEvaporation(G4VE
218 {
219 if(nullptr != ptr && ptr != thePhotonEvapora
220 delete thePhotonEvaporation;
221 thePhotonEvaporation = ptr;
222 if(nullptr != theEvaporation) {
223 theEvaporation->SetPhotonEvaporation(ptr
224 }
225 if(fVerbose > 1) {
226 G4cout << "G4ExcitationHandler::SetPhoto
227 << " for handler " << this << G4e
228 }
229 }
230 }
231
232 void G4ExcitationHandler::SetDeexChannelsType(
233 {
234 G4Evaporation* evap = static_cast<G4Evaporat
235 if(fVerbose > 1) {
236 G4cout << "G4ExcitationHandler::SetDeexCha
237 << " for " << this << G4endl;
238 }
239 if(val == fDummy) {
240 isActive = false;
241 return;
242 }
243 if (nullptr == evap) { return; }
244 if (val == fEvaporation) {
245 evap->SetDefaultChannel();
246 } else if (val == fCombined) {
247 evap->SetCombinedChannel();
248 } else if (val == fGEM) {
249 evap->SetGEMChannel();
250 } else if (val == fGEMVI) {
251 evap->SetGEMVIChannel();
252 }
253 evap->InitialiseChannels();
254 if (fVerbose > 1) {
255 if (G4Threading::IsMasterThread()) {
256 G4cout << "Number of de-excitation chann
257 << theEvaporation->GetNumberOfChannels(
258 G4cout << " " << this;
259 }
260 G4cout << G4endl;
261 }
262 }
263
264 G4VEvaporation* G4ExcitationHandler::GetEvapor
265 {
266 if (nullptr != theEvaporation) { SetParamete
267 return theEvaporation;
268 }
269
270 G4VMultiFragmentation* G4ExcitationHandler::Ge
271 {
272 if (nullptr != theMultiFragmentation) { SetP
273 return theMultiFragmentation;
274 }
275
276 G4VFermiBreakUp* G4ExcitationHandler::GetFermi
277 {
278 if (nullptr != theFermiModel) { SetParameter
279 return theFermiModel;
280 }
281
282 G4VEvaporationChannel* G4ExcitationHandler::Ge
283 {
284 if(nullptr != thePhotonEvaporation) { SetPar
285 return thePhotonEvaporation;
286 }
287
288 G4ReactionProductVector *
289 G4ExcitationHandler::BreakItUp(const G4Fragmen
290 {
291 // Variables existing until end of method
292 G4Fragment * theInitialStatePtr = new G4Frag
293 if (fVerbose > 1) {
294 G4cout << "@@@@@@@@@@ Start G4Excitation H
295 G4cout << theInitialState << G4endl;
296 }
297 if (!isInitialised) { Initialise(); }
298
299 // pointer to fragment vector which receives
300 G4FragmentVector * theTempResult = nullptr;
301
302 theResults.clear();
303 theEvapList.clear();
304
305 // Variables to describe the excited configu
306 G4double exEnergy = theInitialState.GetExcit
307 G4int A = theInitialState.GetA_asInt();
308 G4int Z = theInitialState.GetZ_asInt();
309 G4int nL = theInitialState.GetNumberOfLambda
310
311 // too much excitation
312 if (exEnergy > A*maxExcitation && A > 0) {
313 ++fWarnings;
314 if(fWarnings < 0) {
315 G4ExceptionDescription ed;
316 ed << "High excitation Fragment Z= " <<
317 << " Eex/A(MeV)= " << exEnergy/A;
318 G4Exception("G4ExcitationHandler::BreakI
319 }
320 }
321
322 // for hyper-nuclei subtract lambdas from th
323 G4double lambdaF = 0.0;
324 G4LorentzVector lambdaLV = theInitialStatePt
325 if (0 < nL) {
326
327 // is it a stable hyper-nuclei?
328 if(A >= 3 && A <= 5 && nL <= 2) {
329 G4int pdg = 0;
330 if(3 == A && 1 == nL) {
331 pdg = 1010010030;
332 } else if(5 == A && 2 == Z && 1 == nL) {
333 pdg = 1010020050;
334 } else if(4 == A) {
335 if(1 == Z && 1 == nL) {
336 pdg = 1010010040;
337 } else if(2 == Z && 1 == nL) {
338 pdg = 1010020040;
339 } else if(0 == Z && 2 == nL) {
340 pdg = 1020000040;
341 } else if(1 == Z && 2 == nL) {
342 pdg = 1020010040;
343 }
344 }
345 // initial state is one of hyper-nuclei
346 if (0 < pdg) {
347 const G4ParticleDefinition* part = thePartTa
348 if(nullptr != part) {
349 G4ReactionProduct* theNew = new G4Reaction
350 G4ThreeVector dir = G4ThreeVector( 0.0, 0.
351 if ( lambdaLV.vect().mag() > CLHEP::
352 dir = lambdaLV.vect().unit();
353 }
354 G4double mass = part->GetPDGMass();
355 G4double etot = std::max(lambdaLV.e(), mas
356 dir *= std::sqrt((etot - mass)*(etot
357 theNew->SetMomentum(dir);
358 theNew->SetTotalEnergy(etot);
359 theNew->SetFormationTime(theInitialState.G
360 theNew->SetCreatorModelID(theInitialState.
361 G4ReactionProductVector* v = new G4Reactio
362 v->push_back(theNew);
363 return v;
364 }
365 }
366 }
367 G4double mass = theInitialStatePtr->GetGro
368 lambdaF = nL*(fLambdaMass - CLHEP::neutron
369
370 // de-excitation with neutrons instead of
371 theInitialStatePtr->SetZAandMomentum(lambd
372
373 // 4-momentum not used in de-excitation
374 lambdaLV *= lambdaF;
375 } else if (0 > nL) {
376 ++fWarnings;
377 if(fWarnings < 0) {
378 G4ExceptionDescription ed;
379 ed << "Fragment with negative L: Z=" <<
380 << " Eex/A(MeV)= " << exEnergy/A;
381 G4Exception("G4ExcitationHandler::BreakI
382 }
383 }
384
385 // In case A <= 1 the fragment will not perf
386 if (A <= 1 || !isActive) {
387 theResults.push_back( theInitialStatePtr )
388
389 // check if a fragment is stable
390 } else if (exEnergy < minExcitation && nist-
391 theResults.push_back( theInitialStatePtr )
392
393 // JMQ 150909: first step in de-excitation
394 // Fragments after the first step are stor
395 } else {
396 if ((A<maxAForFermiBreakUp && Z<maxZForFer
397 || exEnergy <= minEForMultiFrag*A) {
398 theEvapList.push_back(theInitialStatePtr
399
400 // Statistical Multifragmentation will tak
401 } else {
402 theTempResult = theMultiFragmentation->B
403 if (nullptr == theTempResult) {
404 theEvapList.push_back(theInitialStatePtr);
405 } else {
406 std::size_t nsec = theTempResult->size();
407
408 // no fragmentation
409 if (0 == nsec) {
410 theEvapList.push_back(theInitialStatePtr);
411
412 // secondary are produced - sort out secon
413 } else {
414 G4bool deletePrimary = true;
415 for (auto const & ptr : *theTempResult) {
416 if (ptr == theInitialStatePtr) { deleteP
417 SortSecondaryFragment(ptr);
418 }
419 if (deletePrimary) { delete theInitialStat
420 }
421 delete theTempResult; // end multifragmentat
422 }
423 }
424 }
425 if (fVerbose > 2) {
426 G4cout << "## After first step of handler
427 << " for evap; "
428 << theResults.size() << " results. " << G
429 }
430 // -----------------------------------
431 // FermiBreakUp and De-excitation loop
432 // -----------------------------------
433
434 static const G4int countmax = 1000;
435 std::size_t kk;
436 for (kk=0; kk<theEvapList.size(); ++kk) {
437 G4Fragment* frag = theEvapList[kk];
438 if (fVerbose > 3) {
439 G4cout << "Next evaporate: " << G4endl;
440 G4cout << *frag << G4endl;
441 }
442 if (kk >= countmax) {
443 G4ExceptionDescription ed;
444 ed << "Infinite loop in the de-excitatio
445 << " iterations \n"
446 << " Initial fragment: \n" << theIniti
447 << "\n Current fragment: \n" << *frag;
448 G4Exception("G4ExcitationHandler::BreakI
449 ed,"Stop execution");
450
451 }
452 A = frag->GetA_asInt();
453 Z = frag->GetZ_asInt();
454 results.clear();
455 if (fVerbose > 2) {
456 G4cout << "G4ExcitationHandler# " << kk
457 << " Eex(MeV)= " << frag->GetExci
458 }
459 // Fermi Break-Up
460 if (theFermiModel->IsApplicable(Z, A, frag
461 theFermiModel->BreakFragment(&results, f
462 std::size_t nsec = results.size();
463 if (fVerbose > 2) { G4cout << "FermiBrea
464
465 // FBU takes care to delete input fragme
466 // The secondary may be excited - photo-
467 if (1 < nsec) {
468 for (auto const & res : results) {
469 SortSecondaryFragment(res);
470 }
471 continue;
472 }
473 // evaporation will be applied
474 }
475 // apply Evaporation, residual nucleus is
476 // photon evaporation is possible
477 theEvaporation->BreakFragment(&results, fr
478 if (fVerbose > 3) {
479 G4cout << "Evaporation Nsec= " << result
480 }
481 if (0 == results.size()) {
482 theResults.push_back(frag);
483 } else {
484 SortSecondaryFragment(frag);
485 }
486
487 // Sort out secondary fragments
488 for (auto const & res : results) {
489 if(fVerbose > 4) {
490 G4cout << "Evaporated product #" << *res <<
491 }
492 SortSecondaryFragment(res);
493 } // end of loop on secondary
494 } // end of the loop over theEvapList
495 if (fVerbose > 2) {
496 G4cout << "## After 2nd step of handler "
497 << " was evap; "
498 << theResults.size() << " results. " << G
499 }
500 G4ReactionProductVector * theReactionProduct
501 new G4ReactionProductVector();
502
503 // MAC (24/07/08)
504 // To optimise the storing speed, we reserve
505 // in memory for the vector
506 theReactionProductVector->reserve( theResult
507
508 if (fVerbose > 1) {
509 G4cout << "### ExcitationHandler provides
510 << " evaporated products:" << G4endl;
511 }
512 G4LorentzVector partOfLambdaLV;
513 if ( nL > 0 ) partOfLambdaLV = lambdaLV/(G4d
514 for (auto const & frag : theResults) {
515 G4LorentzVector lv0 = frag->GetMomentum();
516 G4double etot = lv0.e();
517
518 // in the case of dummy de-excitation, exc
519 // into kinetic energy of output ion
520 if (!isActive) {
521 G4double mass = frag->GetGroundStateMass
522 G4double ptot = lv0.vect().mag();
523 G4double fac = (etot <= mass || 0.0 ==
524 : std::sqrt((etot - mass)*(etot + mass))/pto
525 G4LorentzVector lv((frag->GetMomentum())
526 (frag->GetMomentum()).py()*fac,
527 (frag->GetMomentum()).pz()*fac, etot);
528 frag->SetMomentum(lv);
529 }
530 if (fVerbose > 3) {
531 G4cout << *frag;
532 if (frag->NuclearPolarization()) {
533 G4cout << " " << frag->NuclearPolarization(
534 }
535 G4cout << G4endl;
536 }
537
538 G4int fragmentA = frag->GetA_asInt();
539 G4int fragmentZ = frag->GetZ_asInt();
540 G4double eexc = 0.0;
541 const G4ParticleDefinition* theKindOfFragm
542 G4bool isHyperN = false;
543 if (fragmentA == 0) { // photon or e
544 theKindOfFragment = frag->GetParticleDef
545 } else if (fragmentA == 1 && fragmentZ ==
546 theKindOfFragment = theNeutron;
547 } else if (fragmentA == 1 && fragmentZ ==
548 theKindOfFragment = theProton;
549 } else if (fragmentA == 2 && fragmentZ ==
550 theKindOfFragment = theDeuteron;
551 } else if (fragmentA == 3 && fragmentZ ==
552 theKindOfFragment = theTriton;
553 if(0 < nL) {
554 const G4ParticleDefinition* p = thePar
555 if(nullptr != p) {
556 theKindOfFragment = p;
557 isHyperN = true;
558 --nL;
559 }
560 }
561 } else if (fragmentA == 3 && fragmentZ ==
562 theKindOfFragment = theHe3;
563 } else if (fragmentA == 4 && fragmentZ ==
564 theKindOfFragment = theAlpha;
565 if (0 < nL) {
566 const G4ParticleDefinition* p = thePar
567 if(nullptr != p) {
568 theKindOfFragment = p;
569 isHyperN = true;
570 --nL;
571 }
572 }
573 } else {
574
575 // fragment
576 eexc = frag->GetExcitationEnergy();
577 G4int idxf = frag->GetFloatingLevelNumbe
578 if (eexc < minExcitation) {
579 eexc = 0.0;
580 idxf = 0;
581 }
582
583 theKindOfFragment = theTableOfIons->GetI
584
585 if (fVerbose > 3) {
586 G4cout << "### EXCH: Find ion Z= " << fragme
587 << " A= " << fragmentA
588 << " Eexc(MeV)= " << eexc/MeV << " id
589 << " " << theKindOfFragment->GetParti
590 << G4endl;
591 }
592 }
593 // fragment identified
594 if (nullptr != theKindOfFragment) {
595 G4ReactionProduct * theNew = new G4React
596 if (isHyperN) {
597 G4LorentzVector lv = lv0 + partOfLambd
598 G4ThreeVector dir = lv.vect().unit();
599 G4double mass = theKindOfFragment->Get
600 etot = std::max(lv.e(), mass);
601 G4double ptot = std::sqrt((etot - mass
602 dir *= ptot;
603 theNew->SetMomentum(dir);
604 // remaining not compensated 4-momentum
605 lambdaLV += (lv0 - G4LorentzVector(dir
606 } else {
607 theNew->SetMomentum(lv0.vect());
608 }
609 theNew->SetTotalEnergy(etot);
610 theNew->SetFormationTime(frag->GetCreati
611 theNew->SetCreatorModelID(frag->GetCreat
612 theReactionProductVector->push_back(theN
613
614 // fragment not found out ground state i
615 } else {
616 theKindOfFragment =
617 theTableOfIons->GetIon(fragmentZ,fragmentA,0
618 if (theKindOfFragment) {
619 G4ThreeVector mom(0.0,0.0,0.0);
620 G4double ionmass = theKindOfFragment->GetPDG
621 if (etot <= ionmass) {
622 etot = ionmass;
623 } else {
624 G4double ptot = std::sqrt((etot - ionmass)
625 mom = (frag->GetMomentum().vect().unit())*
626 }
627 G4ReactionProduct * theNew = new G4ReactionP
628 theNew->SetMomentum(mom);
629 theNew->SetTotalEnergy(etot);
630 theNew->SetFormationTime(frag->GetCreationTi
631 theNew->SetCreatorModelID(frag->GetCreatorMo
632 theReactionProductVector->push_back(theNew);
633 if (fVerbose > 3) {
634 G4cout << " ground state, energy
635 << etot << G4endl;
636 }
637 }
638 }
639 delete frag;
640 }
641 // remaining lambdas are free; conserve quan
642 // not 4-momentum
643 if (0 < nL) {
644 G4ThreeVector dir = G4ThreeVector(0.0, 0.0
645 if (lambdaLV.vect().mag() > CLHEP::eV) {
646 dir = lambdaLV.vect().unit();
647 }
648 G4double etot = std::max(lambdaLV.e()/(G4d
649 dir *= std::sqrt((etot - fLambdaMass)*(eto
650 for (G4int i=0; i<nL; ++i) {
651 G4ReactionProduct* theNew = new G4Reacti
652 theNew->SetMomentum(dir);
653 theNew->SetTotalEnergy(etot);
654 theNew->SetFormationTime(theInitialState
655 theNew->SetCreatorModelID(theInitialStat
656 theReactionProductVector->push_back(theN
657 }
658 }
659 if (fVerbose > 3) {
660 G4cout << "@@@@@@@@@@ End G4Excitation Han
661 }
662 return theReactionProductVector;
663 }
664
665 void G4ExcitationHandler::ModelDescription(std
666 {
667 outFile << "G4ExcitationHandler description\
668 << "This class samples de-excitation of ex
669 << "Fermi Break-up model for light fragmen
670 << "evaporation, fission, and photo-evapor
671 << "particle may be proton, neutron, and o
672 << "(Z < 13, A < 29). During photon evapor
673 << "or electrons due to internal conversio
674 }
675
676
677
678