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