Geant4 Cross Reference |
1 // 2 // ******************************************************************** 3 // * License and Disclaimer * 4 // * * 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. * 10 // * * 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitation of liability. * 17 // * * 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************************************** 25 // 26 // 27 // FermiBreakUp de-excitation model 28 // by V. Ivanchenko (July 2016) 29 // 30 31 #include "G4FermiFragmentsPoolVI.hh" 32 #include "G4PhysicalConstants.hh" 33 #include "G4SystemOfUnits.hh" 34 #include "G4NuclearLevelData.hh" 35 #include "G4NucleiProperties.hh" 36 #include "G4LevelManager.hh" 37 #include "G4DeexPrecoParameters.hh" 38 #include "G4FermiBreakUpUtil.hh" 39 #include <iomanip> 40 41 G4FermiFragmentsPoolVI::G4FermiFragmentsPoolVI() 42 {} 43 44 G4FermiFragmentsPoolVI::~G4FermiFragmentsPoolVI() 45 { 46 for (G4int i=0; i<maxA; ++i) { 47 for (G4int j=0; j<maxZ; ++j) { 48 auto ptr = list_c[j][i]; 49 if (nullptr != ptr) { 50 for ( auto const & p : *ptr) { delete p; } 51 delete ptr; 52 } 53 } 54 } 55 for (auto const & ptr : fragment_pool) { delete ptr; } 56 } 57 58 G4bool G4FermiFragmentsPoolVI::HasDecay(const G4int Z, const G4int A, 59 const G4double eexc) const 60 { 61 if (Z < maxZ && A < maxA && nullptr != list_c[Z][A]) { 62 for (auto const & ch : *(list_c[Z][A])) { 63 if (ch->GetExcitation() <= eexc + fTolerance && ch->NumberPairs() > 0) { 64 return true; 65 } 66 } 67 } 68 return false; 69 } 70 71 const G4FermiChannels* 72 G4FermiFragmentsPoolVI::ClosestChannels(G4int Z, G4int A, G4double etot) const 73 { 74 const G4FermiChannels* res = nullptr; 75 if (Z >= maxZ || A >= maxA) { return res; } 76 77 auto chan = list_c[Z][A]; 78 if (nullptr == chan) { return res; } 79 80 G4double demax = 1.e+9; 81 for (auto const & ch : *chan) { 82 if (ch->NumberPairs() == 0) { continue; } 83 G4double de = etot - ch->GetFragment()->GetTotalEnergy(); 84 // an excitation coincide with a level 85 if (std::abs(de) <= fTolerance) { 86 return ch; 87 } 88 if (de >= 0 && de < demax) { 89 demax = de; 90 res = ch; 91 } 92 } 93 return res; 94 } 95 96 G4bool G4FermiFragmentsPoolVI::IsInThePool(const G4int Z, const G4int A, 97 const G4double exc) const 98 { 99 for (auto const& fr : fragment_pool) { 100 if(fr->GetZ() == Z && fr->GetA() == A && 101 std::abs(exc - fr->GetExcitationEnergy()) < fTolerance) 102 { return true; } 103 } 104 return false; 105 } 106 107 void G4FermiFragmentsPoolVI::Initialise() 108 { 109 if (isInitialized) { return; } 110 isInitialized = true; 111 G4DeexPrecoParameters* param = 112 G4NuclearLevelData::GetInstance()->GetParameters(); 113 fTolerance = 2*CLHEP::eV; 114 fElim = param->GetFBUEnergyLimit(); 115 116 fragment_pool.reserve(991); 117 118 // stable particles 119 fragment_pool.push_back(new G4FermiFragment(1, 0, 1, 0.0, DBL_MAX)); 120 fragment_pool.push_back(new G4FermiFragment(1, 1, 1, 0.0, DBL_MAX)); 121 fragment_pool.push_back(new G4FermiFragment(2, 1, 2, 0.0, DBL_MAX)); 122 fragment_pool.push_back(new G4FermiFragment(3, 1, 1, 0.0, 3.8879e+08)); 123 fragment_pool.push_back(new G4FermiFragment(3, 2, 1, 0.0, DBL_MAX)); 124 fragment_pool.push_back(new G4FermiFragment(4, 2, 0, 0.0, DBL_MAX)); 125 fragment_pool.push_back(new G4FermiFragment(5, 2, 3, 0.0, 7.0325e-22)); 126 fragment_pool.push_back(new G4FermiFragment(5, 3, 3, 0.0, 3.70493e-22)); 127 128 // use level data and construct the pool 129 G4NuclearLevelData* ndata = G4NuclearLevelData::GetInstance(); 130 for (G4int Z=1; Z<maxZ; ++Z) { 131 G4int Amin = ndata->GetMinA(Z); 132 G4int Amax = std::min(maxA, ndata->GetMaxA(Z)+1); 133 for (G4int A=Amin; A<Amax; ++A) { 134 const G4LevelManager* man = ndata->GetLevelManager(Z, A); 135 if (nullptr != man) { 136 std::size_t nn = man->NumberOfTransitions(); 137 for(std::size_t i=0; i<=nn; ++i) { 138 G4double exc = man->LevelEnergy(i); 139 /* 140 G4cout << "++ Z=" << Z << " A=" << A << " Eex=" << exc 141 << " time(ns)=" << man->LifeTime(i)/ns << " i=" << i 142 << " InPool=" << IsInThePool(Z, A, exc) << G4endl; 143 */ 144 // only levels below limit are consided 145 if (exc >= fElim) { continue; } 146 // only new are considered 147 if (IsInThePool(Z, A, exc)) { continue; } 148 fragment_pool.push_back(new G4FermiFragment(A, Z, man->TwoSpinParity(i), 149 exc, man->LifeTime(i))); 150 } 151 } 152 } 153 } 154 G4int nfrag = (G4int)fragment_pool.size(); 155 for (auto const& f : fragment_pool) { 156 G4int Z = f->GetZ(); 157 G4int A = f->GetA(); 158 if (list_c[Z][A] == nullptr) { 159 list_c[Z][A] = new std::vector<G4FermiChannels*>; 160 } 161 (list_c[Z][A])->push_back(new G4FermiChannels(f)); 162 } 163 164 // list of fragment pairs ordered by A 165 for (G4int i=0; i<nfrag; ++i) { 166 const G4FermiFragment* f1 = fragment_pool[i]; 167 G4int Z1 = f1->GetZ(); 168 G4int A1 = f1->GetA(); 169 G4double e1 = f1->GetTotalEnergy(); 170 for (G4int j=0; j<nfrag; ++j) { 171 const G4FermiFragment* f2 = fragment_pool[j]; 172 G4int Z2 = f2->GetZ(); 173 G4int A2 = f2->GetA(); 174 if(A2 < A1 || (A2 == A1 && Z2 < Z1)) { continue; } 175 G4int Z = Z1 + Z2; 176 G4int A = A1 + A2; 177 178 if(Z >= maxZ || A >= maxA) { continue; } 179 180 G4double e2 = f2->GetTotalEnergy(); 181 G4double minE = e1 + e2; 182 G4double exc = minE - G4NucleiProperties::GetNuclearMass(A, Z); 183 /* 184 if(1 == Z) { 185 G4cout << "+!+ Z=" << Z << " A=" << A 186 << " Z1=" << Z1 << " A1=" << A1 187 << " Z2=" << Z2 << " A2=" << A2 << " Eex=" << exc 188 << " Qb=" << G4FermiBreakUpUtil::CoulombBarrier(Z1, A1, Z2, A2, exc) 189 << " e1=" << e1 << " e2=" << e2 190 << " M=" << G4NucleiProperties::GetNuclearMass(A, Z) 191 << G4endl; 192 } 193 */ 194 // ignore very excited case 195 if (exc > fElim) { continue; } 196 auto chan = list_c[Z][A]; 197 if (nullptr == chan) { continue; } 198 std::size_t kmax = chan->size(); 199 for (std::size_t k=0; k<kmax; ++k) { 200 auto ch = (*chan)[k]; 201 const G4double e0 = ch->GetMass(); 202 auto f0 = ch->GetFragment(); 203 if (e0 > minE && G4FermiBreakUpUtil::CheckSpinParity(f1, f2, f0)) { 204 const G4double cb = 205 G4FermiBreakUpUtil::CoulombBarrier(Z1, A1, Z2, A2, ch->GetExcitation()); 206 if (e0 >= minE + cb) { 207 ch->AddChannel(new G4FermiPair(f1, f2)); 208 } 209 } 210 } 211 } 212 } 213 // compute cumulative probabilities 214 for (G4int A=1; A<maxA; ++A) { 215 for (G4int Z=0; Z<maxZ; ++Z) { 216 auto chan = list_c[Z][A]; 217 if(nullptr == chan) { continue; } 218 std::size_t kmax = chan->size(); 219 for (std::size_t k=0; k<kmax; ++k) { 220 auto ch = (*chan)[k]; 221 auto frag = ch->GetFragment(); 222 std::size_t nch = ch->NumberPairs(); 223 if (1 < nch) { 224 const std::vector<G4FermiPair*>& pairs = ch->GetChannels(); 225 G4double ptot = 0.0; 226 for (std::size_t i=0; i<nch; ++i) { 227 ptot += G4FermiBreakUpUtil::Probability(frag->GetA(), 228 pairs[i]->GetFragment1(), 229 pairs[i]->GetFragment2(), 230 frag->GetTotalEnergy(), 231 frag->GetExcitationEnergy()); 232 pairs[i]->SetProbability(ptot); 233 } 234 // normalisation 235 if (0.0 == ptot) { 236 pairs[0]->SetProbability(1.0); 237 } else { 238 ptot = 1./ptot; 239 for (std::size_t i=0; i<nch-1; ++i) { 240 G4double x = ptot*pairs[i]->Probability(); 241 pairs[i]->SetProbability(x); 242 } 243 pairs[nch - 1]->SetProbability(1.0); 244 } 245 } 246 } 247 } 248 } 249 } 250 251 void G4FermiFragmentsPoolVI::DumpFragment(const G4FermiFragment* f) const 252 { 253 if (nullptr != f) { 254 G4long prec = G4cout.precision(6); 255 G4cout << " Z=" << f->GetZ() << " A=" << std::setw(2) << f->GetA() 256 << " Mass(GeV)=" << std::setw(8) << f->GetFragmentMass()/GeV 257 << " Eexc(MeV)=" << std::setw(7) << f->GetExcitationEnergy() 258 << " 2S=" << f->TwoSpinParity() << G4endl; 259 G4cout.precision(prec); 260 } 261 } 262 263 void G4FermiFragmentsPoolVI::Dump() const 264 { 265 G4cout <<"----------------------------------------------------------------" 266 <<G4endl; 267 G4cout << "##### List of Fragments in the Fermi Fragment Pool #####" 268 << G4endl; 269 std::size_t nfrag = fragment_pool.size(); 270 G4cout << " Nfragnents=" << nfrag << " Elim(MeV)=" << fElim/CLHEP::MeV << G4endl; 271 for(std::size_t i=0; i<nfrag; ++i) { 272 DumpFragment(fragment_pool[i]); 273 } 274 G4cout << G4endl; 275 276 277 G4cout << "----------------------------------------------------------------" 278 << G4endl; 279 G4cout << "### G4FermiFragmentPoolVI: fragments sorted by A" << G4endl; 280 281 G4int ama{0}; 282 G4long prec = G4cout.precision(6); 283 for (G4int A=1; A<maxA; ++A) { 284 for (G4int Z=0; Z<maxZ; ++Z) { 285 auto chan = list_c[Z][A]; 286 if (nullptr == chan) { continue; } 287 std::size_t jmax = chan->size(); 288 G4cout << " # A=" << A << " Z=" << Z << " Nfagments=" << jmax << G4endl; 289 for(std::size_t j=0; j<jmax; ++j) { 290 auto ch = (*chan)[j]; 291 if(nullptr == ch) { continue; } 292 auto f = ch->GetFragment(); 293 G4int a1 = f->GetA(); 294 G4int z1 = f->GetZ(); 295 std::size_t nch = ch->NumberPairs(); 296 ama += nch; 297 G4cout << " ("<<a1<<","<<z1<<"); Eex(MeV)= " 298 << f->GetExcitationEnergy() 299 << " 2S=" << f->TwoSpinParity() 300 << "; Nchannels=" << nch 301 << G4endl; 302 for (std::size_t k=0; k<nch; ++k) { 303 auto fpair = ch->GetPair(k); 304 if(nullptr == fpair) { continue; } 305 G4cout << " (" << fpair->GetFragment1()->GetZ() 306 << ", " << fpair->GetFragment1()->GetA() 307 << ", " << fpair->GetFragment1()->GetExcitationEnergy() 308 << ") ("<< fpair->GetFragment2()->GetZ() 309 << ", " << fpair->GetFragment2()->GetA() << ", " 310 << fpair->GetFragment2()->GetExcitationEnergy() 311 << ") prob= " << fpair->Probability() 312 << G4endl; 313 } 314 } 315 } 316 } 317 G4cout.precision(prec); 318 G4cout << " ======== Total number of channels " << ama << " ======" << G4endl; 319 } 320 321