304{
305
306 G4Fragment * theInitialStatePtr = new G4Fragment(theInitialState);
307 if (fVerbose > 1) {
308 G4cout <<
"@@@@@@@@@@ Start G4Excitation Handler @@@@@@@@@@@@@ " <<
G4endl;
310 }
312
313 theResults.clear();
314 theEvapList.clear();
315
316
321
322
325 if (0 < nL) {
326
327
328 if(
A >= 3 &&
A <= 5 && nL <= 2) {
330 if(3 ==
A && 1 == nL) {
331 pdg = 1010010030;
332 }
else if(5 ==
A && 2 == Z && 1 == nL) {
333 pdg = 1010020050;
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
346 if (0 < pdg) {
347 const G4ParticleDefinition* part = thePartTable->FindParticle(pdg);
348 if(nullptr != part) {
349 G4ReactionProduct* theNew = new G4ReactionProduct(part);
351 if ( lambdaLV.
vect().
mag() > CLHEP::eV ) {
353 }
355 G4double etot = std::max(lambdaLV.
e(), mass);
356 dir *= std::sqrt((etot - mass)*(etot + mass));
362 v->push_back(theNew);
363 return v;
364 }
365 }
366 }
368 lambdaF = nL*(fLambdaMass - CLHEP::neutron_mass_c2)/mass;
369
370
372
373
374 lambdaLV *= lambdaF;
375 } else if (0 > nL) {
376 ++fWarnings;
377 if(fWarnings < 0) {
379 ed <<
"Fragment with negative L: Z=" << Z <<
" A=" <<
A <<
" L=" << nL
380 <<
" Eex/A(MeV)= " << exEnergy/
A;
382 }
383 }
384
385
387 theResults.push_back( theInitialStatePtr );
388
389
390 }
else if (exEnergy < minExcitation && nist->GetIsotopeAbundance(Z,
A) > 0.0) {
391 theResults.push_back( theInitialStatePtr );
392
393
394
395 } else {
396 theEvapList.push_back(theInitialStatePtr);
397 }
398 if (fVerbose > 2) {
399 G4cout <<
"## After first step of handler " << theEvapList.size()
400 << " for evap; "
401 << theResults.size() <<
" results. " <<
G4endl;
402 }
403
404
405
406
407 static const G4int countmax = 1000;
408 std::size_t kk;
409 for (kk=0; kk<theEvapList.size(); ++kk) {
410 G4Fragment* frag = theEvapList[kk];
411 if (fVerbose > 3) {
414 }
415 if (kk >= countmax) {
417 ed << "Infinite loop in the de-excitation module: " << kk
418 << " iterations \n"
419 << " Initial fragment: \n" << theInitialState
420 << "\n Current fragment: \n" << *frag;
422 ed,"Stop execution");
423
424 }
427 results.clear();
428 if (fVerbose > 2) {
429 G4cout <<
"G4ExcitationHandler# " << kk <<
" Z= " << Z <<
" A= " <<
A
431 }
432
434 theFermiModel->BreakFragment(&results, frag);
435 std::size_t nsec = results.size();
436 if (fVerbose > 2) {
G4cout <<
"FermiBreakUp Nsec= " << nsec <<
G4endl; }
437
438
439
440 if (1 < nsec) {
441 for (auto const & res : results) {
442 SortSecondaryFragment(res);
443 }
444 continue;
445 }
446
447 }
448
449
450 theEvaporation->BreakFragment(&results, frag);
451 if (fVerbose > 3) {
452 G4cout << kk <<
". Evaporation: Nsec=" << results.size()
458 }
459 if (0 == results.size()) {
460 theResults.push_back(frag);
461 } else {
462 SortSecondaryFragment(frag);
463 }
464
465
466 for (auto const & res : results) {
467 if(fVerbose > 4) {
469 }
470 SortSecondaryFragment(res);
471 }
472 }
473 if (fVerbose > 2) {
474 G4cout <<
"## After 2nd step of handler " << theEvapList.size()
475 << " was evap; "
476 << theResults.size() <<
" results. " <<
G4endl;
477 }
480
481
482
483 theReactionProductVector->reserve(theResults.size());
484
485 if (fVerbose > 1) {
486 G4cout <<
"### ExcitationHandler provides " << theResults.size()
487 <<
" evaporated products:" <<
G4endl;
488 }
490 if ( nL > 0 ) partOfLambdaLV = lambdaLV/(
G4double)nL;
491 for (auto const & frag : theResults) {
494
495
496
497 if (!isActive) {
500 G4double fac = (etot <= mass || 0.0 == ptot) ? 0.0
501 : std::sqrt((etot - mass)*(etot + mass))/ptot;
506 }
507 if (fVerbose > 3) {
509 if (frag->NuclearPolarization()) {
510 G4cout <<
" " << frag->NuclearPolarization();
511 }
513 }
514
515 G4int fragmentA = frag->GetA_asInt();
516 G4int fragmentZ = frag->GetZ_asInt();
518 const G4ParticleDefinition* theKindOfFragment = nullptr;
520 if (fragmentA == 0) {
521 theKindOfFragment = frag->GetParticleDefinition();
522 } else if (fragmentA == 1 && fragmentZ == 0) {
523 theKindOfFragment = theNeutron;
524 } else if (fragmentA == 1 && fragmentZ == 1) {
525 theKindOfFragment = theProton;
526 } else if (fragmentA == 2 && fragmentZ == 1) {
527 theKindOfFragment = theDeuteron;
528 } else if (fragmentA == 3 && fragmentZ == 1) {
529 theKindOfFragment = theTriton;
530 if(0 < nL) {
531 const G4ParticleDefinition* p = thePartTable->FindParticle(1010010030);
532 if(nullptr != p) {
533 theKindOfFragment = p;
534 isHyperN = true;
535 --nL;
536 }
537 }
538 } else if (fragmentA == 3 && fragmentZ == 2) {
539 theKindOfFragment = theHe3;
540 } else if (fragmentA == 4 && fragmentZ == 2) {
541 theKindOfFragment = theAlpha;
542 if (0 < nL) {
543 const G4ParticleDefinition* p = thePartTable->FindParticle(1010020040);
544 if(nullptr != p) {
545 theKindOfFragment = p;
546 isHyperN = true;
547 --nL;
548 }
549 }
550 } else {
551
552
553 eexc = frag->GetExcitationEnergy();
554 G4int idxf = frag->GetFloatingLevelNumber();
555 if (eexc < minExcitation) {
556 eexc = 0.0;
557 idxf = 0;
558 }
559
560 theKindOfFragment = theTableOfIons->GetIon(fragmentZ, fragmentA, eexc,
562 if (fVerbose > 3) {
563 G4cout <<
"### EXCH: Find ion Z= " << fragmentZ
564 << " A= " << fragmentA
565 << " Eexc(MeV)= " << eexc/MeV << " idx= " << idxf
568 }
569 }
570
571 if (nullptr != theKindOfFragment) {
572 G4ReactionProduct * theNew = new G4ReactionProduct(theKindOfFragment);
573 if (isHyperN) {
577 etot = std::max(lv.
e(), mass);
578 G4double ptot = std::sqrt((etot - mass)*(etot + mass));
579 dir *= ptot;
581
583 } else {
585 }
589 theReactionProductVector->push_back(theNew);
590
591
592 } else {
593 theKindOfFragment =
594 theTableOfIons->GetIon(fragmentZ,fragmentA,0.0,
noFloat,0);
595 if (theKindOfFragment) {
598 if (etot <= ionmass) {
599 etot = ionmass;
600 } else {
601 G4double ptot = std::sqrt((etot - ionmass)*(etot + ionmass));
602 mom = (frag->GetMomentum().vect().unit())*ptot;
603 }
604 G4ReactionProduct * theNew = new G4ReactionProduct(theKindOfFragment);
609 theReactionProductVector->push_back(theNew);
610 if (fVerbose > 3) {
611 G4cout <<
" ground state, energy corrected E(MeV)= "
613 }
614 }
615 }
616 delete frag;
617 }
618
619
620 if (0 < nL) {
622 if (lambdaLV.
vect().
mag() > CLHEP::eV) {
624 }
626 dir *= std::sqrt((etot - fLambdaMass)*(etot + fLambdaMass));
627 for (
G4int i=0; i<nL; ++i) {
628 G4ReactionProduct* theNew = new G4ReactionProduct(theLambda);
633 theReactionProductVector->push_back(theNew);
634 }
635 }
636 if (fVerbose > 3) {
637 G4cout <<
"@@@@@@@@@@ End G4Excitation Handler "<<
G4endl;
638 }
639 return theReactionProductVector;
640}
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
std::ostringstream G4ExceptionDescription
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
G4bool IsLongLived() const
G4int GetNumberOfLambdas() const
void SetMomentum(const G4LorentzVector &value)
void SetZAandMomentum(const G4LorentzVector &, G4int Z, G4int A, G4int nLambdas=0)
static G4Ions::G4FloatLevelBase FloatLevelBase(char flbChar)
G4double GetPDGMass() const
const G4String & GetParticleName() const
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)