291{
292
294 if (fVerbose > 1) {
295 G4cout <<
"@@@@@@@@@@ Start G4Excitation Handler @@@@@@@@@@@@@ " <<
G4endl;
297 }
299
300
302
303 theResults.clear();
304 theEvapList.clear();
305
306
311
312
313 if (exEnergy >
A*maxExcitation &&
A > 0) {
314 ++fWarnings;
315 if(fWarnings < 0) {
317 ed <<
"High excitation Fragment Z= " << Z <<
" A= " <<
A
318 <<
" Eex/A(MeV)= " << exEnergy/
A;
320 }
321 }
322
323
326 if (0 < nL) {
327
328
329 if(
A >= 3 &&
A <= 5 && nL <= 2) {
331 if(3 ==
A && 1 == nL) {
332 pdg = 1010010030;
333 }
else if(5 ==
A && 2 == Z && 1 == nL) {
334 pdg = 1010020050;
336 if(1 == Z && 1 == nL) {
337 pdg = 1010010040;
338 } else if(2 == Z && 1 == nL) {
339 pdg = 1010020040;
340 } else if(0 == Z && 2 == nL) {
341 pdg = 1020000040;
342 } else if(1 == Z && 2 == nL) {
343 pdg = 1020010040;
344 }
345 }
346
347 if (0 < pdg) {
349 if(nullptr != part) {
352 if ( lambdaLV.
vect().
mag() > CLHEP::eV ) {
354 }
356 G4double etot = std::max(lambdaLV.
e(), mass);
357 dir *= std::sqrt((etot - mass)*(etot + mass));
363 v->push_back(theNew);
364 return v;
365 }
366 }
367 }
369 lambdaF = nL*(fLambdaMass - CLHEP::neutron_mass_c2)/mass;
370
371
373
374
375 lambdaLV *= lambdaF;
376 } else if (0 > nL) {
377 ++fWarnings;
378 if(fWarnings < 0) {
380 ed <<
"Fragment with negative L: Z=" << Z <<
" A=" <<
A <<
" L=" << nL
381 <<
" Eex/A(MeV)= " << exEnergy/
A;
383 }
384 }
385
386
387 if (
A <= 1 || !isActive) {
388 theResults.push_back( theInitialStatePtr );
389
390
391 }
else if (exEnergy < minExcitation && nist->GetIsotopeAbundance(Z,
A) > 0.0) {
392 theResults.push_back( theInitialStatePtr );
393
394
395
396 } else {
397 if ((
A<maxAForFermiBreakUp && Z<maxZForFermiBreakUp)
398 || exEnergy <= minEForMultiFrag*
A) {
399 theEvapList.push_back(theInitialStatePtr);
400
401
402 } else {
403 theTempResult = theMultiFragmentation->
BreakItUp(theInitialState);
404 if (nullptr == theTempResult) {
405 theEvapList.push_back(theInitialStatePtr);
406 } else {
407 std::size_t nsec = theTempResult->size();
408
409
410 if (0 == nsec) {
411 theEvapList.push_back(theInitialStatePtr);
412
413
414 } else {
415 G4bool deletePrimary =
true;
416 for (auto const & ptr : *theTempResult) {
417 if (ptr == theInitialStatePtr) { deletePrimary = false; }
418 SortSecondaryFragment(ptr);
419 }
420 if (deletePrimary) { delete theInitialStatePtr; }
421 }
422 delete theTempResult;
423 }
424 }
425 }
426 if (fVerbose > 2) {
427 G4cout <<
"## After first step of handler " << theEvapList.size()
428 << " for evap; "
429 << theResults.size() <<
" results. " <<
G4endl;
430 }
431
432
433
434
435 static const G4int countmax = 1000;
436 std::size_t kk;
437 for (kk=0; kk<theEvapList.size(); ++kk) {
439 if (fVerbose > 3) {
442 }
443 if (kk >= countmax) {
445 ed << "Infinite loop in the de-excitation module: " << kk
446 << " iterations \n"
447 << " Initial fragment: \n" << theInitialState
448 << "\n Current fragment: \n" << *frag;
450 ed,"Stop execution");
451
452 }
455 results.clear();
456 if (fVerbose > 2) {
457 G4cout <<
"G4ExcitationHandler# " << kk <<
" Z= " << Z <<
" A= " <<
A
459 }
460
463 std::size_t nsec = results.size();
464 if (fVerbose > 2) {
G4cout <<
"FermiBreakUp Nsec= " << nsec <<
G4endl; }
465
466
467
468 if (1 < nsec) {
469 for (auto const & res : results) {
470 SortSecondaryFragment(res);
471 }
472 continue;
473 }
474
475 }
476
477
479 if (fVerbose > 3) {
480 G4cout <<
"Evaporation Nsec= " << results.size() <<
G4endl;
481 }
482 if (0 == results.size()) {
483 theResults.push_back(frag);
484 } else {
485 SortSecondaryFragment(frag);
486 }
487
488
489 for (auto const & res : results) {
490 if(fVerbose > 4) {
492 }
493 SortSecondaryFragment(res);
494 }
495 }
496 if (fVerbose > 2) {
497 G4cout <<
"## After 2nd step of handler " << theEvapList.size()
498 << " was evap; "
499 << theResults.size() <<
" results. " <<
G4endl;
500 }
503
504
505
506
507 theReactionProductVector->reserve( theResults.size() );
508
509 if (fVerbose > 2) {
510 G4cout <<
"### ExcitationHandler provides " << theResults.size()
511 <<
" evaporated products:" <<
G4endl;
512 }
514 if ( nL > 0 ) partOfLambdaLV = lambdaLV/(
G4double)nL;
515 for (auto const & frag : theResults) {
518
519
520
521 if (!isActive) {
524 G4double fac = (etot <= mass || 0.0 == ptot) ? 0.0
525 :
std::sqrt((etot - mass)*(etot + mass))/ptot;
530 }
531 if (fVerbose > 3) {
533 if (frag->NuclearPolarization()) {
534 G4cout <<
" " << frag->NuclearPolarization();
535 }
537 }
538
539 G4int fragmentA = frag->GetA_asInt();
540 G4int fragmentZ = frag->GetZ_asInt();
544 if (fragmentA == 0) {
545 theKindOfFragment = frag->GetParticleDefinition();
546 } else if (fragmentA == 1 && fragmentZ == 0) {
547 theKindOfFragment = theNeutron;
548 } else if (fragmentA == 1 && fragmentZ == 1) {
549 theKindOfFragment = theProton;
550 } else if (fragmentA == 2 && fragmentZ == 1) {
551 theKindOfFragment = theDeuteron;
552 } else if (fragmentA == 3 && fragmentZ == 1) {
553 theKindOfFragment = theTriton;
554 if(0 < nL) {
556 if(nullptr != p) {
557 theKindOfFragment = p;
558 isHyperN = true;
559 --nL;
560 }
561 }
562 } else if (fragmentA == 3 && fragmentZ == 2) {
563 theKindOfFragment = theHe3;
564 } else if (fragmentA == 4 && fragmentZ == 2) {
565 theKindOfFragment = theAlpha;
566 if (0 < nL) {
568 if(nullptr != p) {
569 theKindOfFragment = p;
570 isHyperN = true;
571 --nL;
572 }
573 }
574 } else {
575
576
577 eexc = frag->GetExcitationEnergy();
578 G4int idxf = frag->GetFloatingLevelNumber();
579 if (eexc < minExcitation) {
580 eexc = 0.0;
581 idxf = 0;
582 }
583
584 theKindOfFragment = theTableOfIons->
GetIon(fragmentZ, fragmentA, eexc,
586 if (fVerbose > 3) {
587 G4cout <<
"### EXCH: Find ion Z= " << fragmentZ
588 << " A= " << fragmentA
589 << " Eexc(MeV)= " << eexc/MeV << " idx= " << idxf
591 }
592 }
593
594 if (nullptr != theKindOfFragment) {
596 if (isHyperN) {
600 etot = std::max(lv.
e(), mass);
601 G4double ptot = std::sqrt((etot - mass)*(etot + mass));
602 dir *= ptot;
604
606 } else {
608 }
611 if (theKindOfFragment == theElectron) {
613 } else {
615 }
616 theReactionProductVector->push_back(theNew);
617
618
619 } else {
620 theKindOfFragment =
622 if (theKindOfFragment) {
625 if (etot <= ionmass) {
626 etot = ionmass;
627 } else {
628 G4double ptot = std::sqrt((etot - ionmass)*(etot + ionmass));
629 mom = (frag->GetMomentum().vect().unit())*ptot;
630 }
636 theReactionProductVector->push_back(theNew);
637 if (fVerbose > 3) {
638 G4cout <<
" ground state, energy corrected E(MeV)= "
640 }
641 }
642 }
643 delete frag;
644 }
645
646
647 if (0 < nL) {
649 if (lambdaLV.
vect().
mag() > CLHEP::eV) {
651 }
653 dir *= std::sqrt((etot - fLambdaMass)*(etot + fLambdaMass));
654 for (
G4int i=0; i<nL; ++i) {
660 theReactionProductVector->push_back(theNew);
661 }
662 }
663 if (fVerbose > 3) {
664 G4cout <<
"@@@@@@@@@@ End G4Excitation Handler "<<
G4endl;
665 }
666 return theReactionProductVector;
667}
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
std::ostringstream G4ExceptionDescription
std::vector< G4Fragment * > G4FragmentVector
CLHEP::HepLorentzVector G4LorentzVector
std::vector< G4ReactionProduct * > G4ReactionProductVector
CLHEP::Hep3Vector G4ThreeVector
G4double GetGroundStateMass() const
G4int GetCreatorModelID() const
G4double GetExcitationEnergy() const
const G4LorentzVector & GetMomentum() const
G4double GetCreationTime() const
G4int GetNumberOfLambdas() const
void SetMomentum(const G4LorentzVector &value)
void SetZAandMomentum(const G4LorentzVector &, G4int Z, G4int A, G4int nLambdas=0)
G4ParticleDefinition * GetIon(G4int Z, G4int A, G4int lvl=0)
static G4Ions::G4FloatLevelBase FloatLevelBase(char flbChar)
G4ParticleDefinition * FindParticle(G4int PDGEncoding)
void SetMomentum(const G4double x, const G4double y, const G4double z)
void SetTotalEnergy(const G4double en)
void SetCreatorModelID(const G4int mod)
void SetFormationTime(G4double aTime)
virtual void BreakFragment(G4FragmentVector *, G4Fragment *theNucleus)
virtual G4bool IsApplicable(G4int Z, G4int A, G4double eexc) const =0
virtual void BreakFragment(G4FragmentVector *results, G4Fragment *theNucleus)=0
virtual G4FragmentVector * BreakItUp(const G4Fragment &theNucleus)=0