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 // Hadronic Process: Nuclear De-excitations 27 // by V. Lara (Oct 1998) 28 // 29 // Modifications: 30 // 28.10.2010 V.Ivanchenko defined members in constructor and cleaned up 31 32 #include "G4VEmissionProbability.hh" 33 #include "G4NuclearLevelData.hh" 34 #include "G4LevelManager.hh" 35 #include "G4DeexPrecoParameters.hh" 36 #include "Randomize.hh" 37 #include "G4Pow.hh" 38 #include "G4Log.hh" 39 #include "G4Exp.hh" 40 41 G4VEmissionProbability::G4VEmissionProbability(G4int Z, G4int A) 42 : pVerbose(1), theZ(Z), theA(A), elimit(CLHEP::MeV) 43 { 44 pNuclearLevelData = G4NuclearLevelData::GetInstance(); 45 pG4pow = G4Pow::GetInstance(); 46 if(A > 0) { pEvapMass = G4NucleiProperties::GetNuclearMass(theA, theZ); } 47 G4DeexPrecoParameters* param = pNuclearLevelData->GetParameters(); 48 OPTxs = param->GetDeexModelType(); 49 } 50 51 void G4VEmissionProbability::Initialise() 52 { 53 G4DeexPrecoParameters* param = pNuclearLevelData->GetParameters(); 54 pVerbose = param->GetVerbose(); 55 fFD = param->GetDiscreteExcitationFlag(); 56 pTolerance = param->GetMinExcitation(); 57 pWidth = param->GetNuclearLevelWidth(); 58 } 59 60 void G4VEmissionProbability::ResetIntegrator(size_t, G4double de, G4double eps) 61 { 62 if(de > 0.0) { elimit = de; } 63 if(eps > 0.0) { accuracy = eps; } 64 } 65 66 G4double G4VEmissionProbability::EmissionProbability(const G4Fragment&, G4double) 67 { 68 return 0.0; 69 } 70 71 G4double G4VEmissionProbability::ComputeProbability(G4double, G4double) 72 { 73 return 0.0; 74 } 75 76 G4double G4VEmissionProbability::IntegrateProbability(G4double elow, 77 G4double ehigh, 78 G4double cb) 79 { 80 pProbability = 0.0; 81 if(elow >= ehigh) { return pProbability; } 82 83 emin = elow; 84 emax = ehigh; 85 eCoulomb = cb; 86 87 const G4double edeltamin = 0.1*CLHEP::MeV; 88 const G4double edeltamax = 2*CLHEP::MeV; 89 G4double edelta = std::min(std::min(elimit, edeltamax), edeltamin); 90 G4double xbin = (emax - emin)/edelta + 1.0; 91 G4int ibin = std::max((G4int)xbin, 4); 92 93 // providing smart binning 94 G4int nbin = ibin*5; 95 edelta = (emax - emin)/ibin; 96 97 G4double x(emin), y(0.0); 98 G4double edelmicro = edelta*0.02; 99 probmax = ComputeProbability(x + edelmicro, eCoulomb); 100 G4double problast = probmax; 101 if(pVerbose > 1) { 102 G4cout << "### G4VEmissionProbability::IntegrateProbability: " 103 << "probmax=" << probmax << " Emin=" << emin 104 << " Emax=" << emax << " QB=" << cb << " nbin=" << nbin 105 << G4endl; 106 } 107 fE1 = fE2 = fP2 = 0.0; 108 G4double emax0 = emax - edelmicro; 109 G4bool endpoint = false; 110 for(G4int i=0; i<nbin; ++i) { 111 x += edelta; 112 if(x >= emax0) { 113 x = emax0; 114 endpoint = true; 115 } 116 y = ComputeProbability(x, eCoulomb); 117 if(pVerbose > 2) { 118 G4cout << " " << i << ". E= " << x << " prob= " << y 119 << " Edel= " << edelta << G4endl; 120 } 121 if(y >= probmax) { 122 probmax = y; 123 } else if(0.0 == fE1 && 2*y < probmax) { 124 fE1 = x; 125 } 126 127 G4double del = (y + problast)*edelta*0.5; 128 pProbability += del; 129 // end of the loop 130 if(del < accuracy*pProbability || endpoint) { break; } 131 problast = y; 132 133 // smart step definition 134 if(del != pProbability && del > 0.8*pProbability && 135 0.7*edelta > edeltamin) { 136 edelta *= 0.7; 137 } else if(del < 0.1*pProbability && 1.5*edelta < edeltamax) { 138 edelta *= 1.5; 139 } 140 } 141 if(fE1 > emin && fE1 < emax) { 142 fE2 = std::max(0.5*(fE1 + emax), emax - edelta); 143 fP2 = 2*ComputeProbability(fE2, eCoulomb); 144 } 145 146 if(pVerbose > 1) { 147 G4cout << " Probability= " << pProbability << " probmax= " 148 << probmax << " emin=" << emin << " emax=" << emax 149 << " E1=" << fE1 << " E2=" << fE2 << G4endl; 150 } 151 return pProbability; 152 } 153 154 G4double G4VEmissionProbability::SampleEnergy() 155 { 156 static const G4double fact = 1.05; 157 static const G4double alim = 0.05; 158 static const G4double blim = 20.; 159 probmax *= fact; 160 161 // two regions with flat and exponential majorant 162 G4double del = emax - emin; 163 G4double p1 = 1.0; 164 G4double p2 = 0.0; 165 G4double a0 = 0.0; 166 G4double a1 = 1.0; 167 G4double x; 168 if(fE1 > 0.0 && fP2 > 0.0 && fP2 < 0.5*probmax) { 169 a0 = G4Log(probmax/fP2)/(fE2 - fE1); 170 del= fE1 - emin; 171 p1 = del; 172 x = a0*(emax - fE1); 173 if(x < blim) { 174 a1 = (x > alim) ? 1.0 - G4Exp(-x) : x*(1.0 - 0.5*x); 175 } 176 p2 = a1/a0; 177 p1 /= (p1 + p2); 178 p2 = 1.0 - p1; 179 } 180 181 if(pVerbose > 1) { 182 G4cout << "### G4VEmissionProbability::SampleEnergy: " 183 << " Emin= " << emin << " Emax= " << emax 184 << "/n E1=" << fE1 << " p1=" << p1 185 << " probmax=" << probmax << " P2=" << fP2 << G4endl; 186 } 187 188 CLHEP::HepRandomEngine* rndm = G4Random::getTheEngine(); 189 const G4int nmax = 1000; 190 G4double ekin, gg, gmax; 191 G4int n = 0; 192 do { 193 ++n; 194 G4double q = rndm->flat(); 195 if (p2 == 0.0) { 196 gmax = probmax; 197 ekin = del*q + emin; 198 } else if (q <= p1) { 199 gmax = probmax; 200 ekin = del*q/p1 + emin; 201 } else { 202 ekin = fE1 - G4Log(1.0 - (q - p1)*a1/p2)/a0; 203 x = a0*(ekin - fE1); 204 gmax = fP2; 205 if(x < blim) { 206 gmax = probmax*((x > alim) ? G4Exp(-x) : 1.0 - x*(1.0 - 0.5*x)); 207 } 208 } 209 gg = ComputeProbability(ekin, eCoulomb); 210 if(pVerbose > 2) { 211 G4cout << " " << n 212 << ". prob= " << gg << " probmax= " << probmax 213 << " Ekin= " << ekin << G4endl; 214 } 215 if((gg > gmax || n > nmax) && pVerbose > 1) { 216 G4cout << "### G4VEmissionProbability::SampleEnergy for Z= " << theZ 217 << " A= " << theA << " Eex(MeV)=" << fExc << " p1=" << p1 218 << "\n Warning n= " << n 219 << " prob/gmax=" << gg/gmax 220 << " prob=" << gg << " gmax=" << gmax << " probmax=" << probmax 221 << "\n Ekin= " << ekin << " Emin= " << emin 222 << " Emax= " << emax << G4endl; 223 } 224 } while(gmax*rndm->flat() > gg && n < nmax); 225 G4double enew = FindRecoilExcitation(ekin); 226 if(pVerbose > 1) { 227 G4cout << "### SampleEnergy: Efinal= " 228 << enew << " E=" << ekin << " Eexc=" << fExcRes << G4endl; 229 } 230 return enew; 231 } 232 233 G4double G4VEmissionProbability::FindRecoilExcitation(const G4double e) 234 { 235 G4double mass = pEvapMass + fExc; 236 237 G4double m02 = pMass*pMass; 238 G4double m12 = mass*mass; 239 G4double m22 = pResMass*pResMass; 240 G4double mres = std::sqrt(m02 + m12 - 2.*pMass*(mass + e)); 241 242 fExcRes = mres - pResMass; 243 244 if(pVerbose > 1) { 245 G4cout << "### FindRecoilExcitation for resZ= " 246 << resZ << " resA= " << resA 247 << " evaporated Z= " << theZ << " A= " << theA 248 << " Ekin= " << e << " Eexc= " << fExcRes << G4endl; 249 } 250 251 // residual nucleus is in the ground state 252 if(fExcRes < pTolerance) { 253 fExcRes = 0.0; 254 return std::max(0.5*(m02 + m12 - m22)/pMass - mass, 0.0); 255 } 256 if(!fFD) { return e; } 257 258 // select final state excitation 259 auto lManager = pNuclearLevelData->GetLevelManager(resZ, resA); 260 if(nullptr == lManager) { return e; } 261 262 // levels are not known 263 if(fExcRes > lManager->MaxLevelEnergy() + pTolerance) { return e; } 264 265 // find level 266 std::size_t idx = lManager->NearestLevelIndex(fExcRes); 267 auto level = lManager->GetLevel(idx); 268 269 // unstable level 270 if (level->GetTimeGamma() == 0.0) { return e; } 271 272 // is possible to use level energy? 273 G4double elevel = lManager->LevelEnergy(idx); 274 if (std::abs(elevel - fExcRes) > pWidth || pMass < mass + pResMass + elevel) { 275 return e; 276 } 277 278 // long-lived level 279 G4double massR = pResMass + elevel; 280 G4double mr2 = massR*massR; 281 fExcRes = elevel; 282 return std::max(0.5*(m02 + m12 - mr2)/pMass - mass, 0.0); 283 } 284