Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************* 24 // ******************************************************************** 25 // 25 // >> 26 // >> 27 // $Id: G4VEmissionProbability.cc 66241 2012-12-13 18:34:42Z gunter $ >> 28 // 26 // Hadronic Process: Nuclear De-excitations 29 // Hadronic Process: Nuclear De-excitations 27 // by V. Lara (Oct 1998) 30 // by V. Lara (Oct 1998) 28 // 31 // 29 // Modifications: 32 // Modifications: 30 // 28.10.2010 V.Ivanchenko defined members in 33 // 28.10.2010 V.Ivanchenko defined members in constructor and cleaned up 31 34 32 #include "G4VEmissionProbability.hh" 35 #include "G4VEmissionProbability.hh" 33 #include "G4NuclearLevelData.hh" 36 #include "G4NuclearLevelData.hh" 34 #include "G4LevelManager.hh" << 35 #include "G4DeexPrecoParameters.hh" 37 #include "G4DeexPrecoParameters.hh" 36 #include "Randomize.hh" 38 #include "Randomize.hh" 37 #include "G4Pow.hh" << 38 #include "G4Log.hh" << 39 #include "G4Exp.hh" << 40 39 41 G4VEmissionProbability::G4VEmissionProbability 40 G4VEmissionProbability::G4VEmissionProbability(G4int Z, G4int A) 42 : pVerbose(1), theZ(Z), theA(A), elimit(CLHE << 41 :OPTxs(3),fVerbose(0),theZ(Z),theA(A), >> 42 LevelDensity(0.1),elimit(CLHEP::MeV),accuracy(0.02) 43 { 43 { 44 pNuclearLevelData = G4NuclearLevelData::GetI << 44 fG4pow = G4Pow::GetInstance(); 45 pG4pow = G4Pow::GetInstance(); << 45 fPairCorr = G4NuclearLevelData::GetInstance()->GetPairingCorrection(); 46 if(A > 0) { pEvapMass = G4NucleiProperties:: << 46 length = nfilled = 0; 47 G4DeexPrecoParameters* param = pNuclearLevel << 47 emin = emax = eCoulomb = probmax = eprobmax = totProbability = 0.0; 48 OPTxs = param->GetDeexModelType(); << 49 } 48 } 50 49 51 void G4VEmissionProbability::Initialise() << 50 G4VEmissionProbability::~G4VEmissionProbability() 52 { << 51 {} 53 G4DeexPrecoParameters* param = pNuclearLevel << 54 pVerbose = param->GetVerbose(); << 55 fFD = param->GetDiscreteExcitationFlag(); << 56 pTolerance = param->GetMinExcitation(); << 57 pWidth = param->GetNuclearLevelWidth(); << 58 } << 59 52 60 void G4VEmissionProbability::ResetIntegrator(s << 53 void G4VEmissionProbability::Initialise() 61 { 54 { 62 if(de > 0.0) { elimit = de; } << 55 G4DeexPrecoParameters* param = G4NuclearLevelData::GetInstance()->GetParameters(); 63 if(eps > 0.0) { accuracy = eps; } << 56 OPTxs = param->GetDeexModelType(); >> 57 LevelDensity = param->GetLevelDensity(); 64 } 58 } 65 59 66 G4double G4VEmissionProbability::EmissionProba << 60 void G4VEmissionProbability::ResetIntegrator(size_t nbin, G4double de, G4double eps) 67 { 61 { 68 return 0.0; << 62 if(nbin > 0) { >> 63 fProb.clear(); >> 64 fEner.clear(); >> 65 fEner.resize(nbin+1, 0.0); >> 66 fProb.resize(nbin+1, 0.0); >> 67 length = nbin; >> 68 } >> 69 if(de > 0.0) { elimit = de; } >> 70 if(accuracy > 0.0) { accuracy = eps; } 69 } 71 } 70 72 71 G4double G4VEmissionProbability::ComputeProbab 73 G4double G4VEmissionProbability::ComputeProbability(G4double, G4double) 72 { 74 { 73 return 0.0; 75 return 0.0; 74 } 76 } 75 77 76 G4double G4VEmissionProbability::IntegrateProb << 78 G4double 77 << 79 G4VEmissionProbability::IntegrateProbability(G4double elow, G4double ehigh, G4double cb) 78 << 79 { 80 { 80 pProbability = 0.0; << 81 totProbability = 0.0; 81 if(elow >= ehigh) { return pProbability; } << 82 if(elow >= ehigh) { return totProbability; } 82 83 83 emin = elow; 84 emin = elow; 84 emax = ehigh; 85 emax = ehigh; 85 eCoulomb = cb; 86 eCoulomb = cb; 86 87 87 const G4double edeltamin = 0.1*CLHEP::MeV; << 88 size_t nbin = (size_t)((emax - emin)/elimit) + 1; 88 const G4double edeltamax = 2*CLHEP::MeV; << 89 G4double edelta = elimit; 89 G4double edelta = std::min(std::min(elimit, << 90 if(nbin < 3) { 90 G4double xbin = (emax - emin)/edelta + 1.0; << 91 nbin = 3; 91 G4int ibin = std::max((G4int)xbin, 4); << 92 edelta = (emax - emin)/(G4double)nbin; 92 << 93 } 93 // providing smart binning << 94 if(nbin > length) { 94 G4int nbin = ibin*5; << 95 fEner.resize(nbin + 1); 95 edelta = (emax - emin)/ibin; << 96 fProb.resize(nbin + 1); 96 << 97 length = nbin; 97 G4double x(emin), y(0.0); << 98 } 98 G4double edelmicro = edelta*0.02; << 99 99 probmax = ComputeProbability(x + edelmicro, << 100 G4double x(emin), del, y; 100 G4double problast = probmax; << 101 G4double problast = ComputeProbability(x, eCoulomb); 101 if(pVerbose > 1) { << 102 eprobmax= emin; 102 G4cout << "### G4VEmissionProbability::Int << 103 probmax = problast; 103 << "probmax=" << probmax << " Emin=" << e << 104 fEner[0] = emin; 104 << " Emax=" << emax << " QB=" << cb << " << 105 fProb[0] = problast; >> 106 size_t i(0); >> 107 if(fVerbose > 1) { >> 108 G4cout << "### G4VEmissionProbability::IntegrateProbability: " >> 109 << " Emax= " << emax << " QB= " << cb << " nbin= " << nbin 105 << G4endl; 110 << G4endl; >> 111 G4cout << " 0. E= " << emin << " prob= " << problast << G4endl; 106 } 112 } 107 fE1 = fE2 = fP2 = 0.0; << 113 static const G4double edeltamin = 0.2*CLHEP::MeV; 108 G4double emax0 = emax - edelmicro; << 114 static const G4double edeltamax = 2*CLHEP::MeV; 109 G4bool endpoint = false; << 115 G4bool peak = true; 110 for(G4int i=0; i<nbin; ++i) { << 116 do { >> 117 ++i; 111 x += edelta; 118 x += edelta; 112 if(x >= emax0) { << 119 if(x > emax) { 113 x = emax0; << 120 x = emax; 114 endpoint = true; << 121 edelta = emax - fEner[i-1]; 115 } 122 } 116 y = ComputeProbability(x, eCoulomb); 123 y = ComputeProbability(x, eCoulomb); 117 if(pVerbose > 2) { << 124 if(fVerbose > 1) { 118 G4cout << " " << i << ". E= " << x < << 125 G4cout << " " << i << ". E= " << x << " prob= " << y 119 << " Edel= " << edelta << G4endl; 126 << " Edel= " << edelta << G4endl; 120 } 127 } 121 if(y >= probmax) { << 128 fEner[i] = x; 122 probmax = y; << 129 fProb[i] = y; 123 } else if(0.0 == fE1 && 2*y < probmax) { << 130 del = (y + problast)*edelta*0.5; 124 fE1 = x; << 131 totProbability += del; 125 } << 126 << 127 G4double del = (y + problast)*edelta*0.5; << 128 pProbability += del; << 129 // end of the loop << 130 if(del < accuracy*pProbability || endpoint << 131 problast = y; << 132 132 133 // smart step definition 133 // smart step definition 134 if(del != pProbability && del > 0.8*pProba << 134 if(del < accuracy*totProbability) { break; } 135 0.7*edelta > edeltamin) { << 135 else if(del != totProbability && del > 0.8*totProbability 136 edelta *= 0.7; << 136 && edelta > edeltamin) 137 } else if(del < 0.1*pProbability && 1.5*ed << 137 { edelta *= 0.5; } 138 edelta *= 1.5; << 138 else if(del < 0.1*totProbability && edelta < edeltamax) 139 } << 139 { edelta *= 2.0; } 140 } << 140 if(y > probmax) { 141 if(fE1 > emin && fE1 < emax) { << 141 probmax = y; 142 fE2 = std::max(0.5*(fE1 + emax), emax - ed << 142 eprobmax = x; 143 fP2 = 2*ComputeProbability(fE2, eCoulomb); << 143 } else if(peak && y < probmax) { 144 } << 144 edelta *= 2.0; >> 145 peak = false; >> 146 } >> 147 problast = y; >> 148 // Loop checking, 10-Mar-2017, Vladimir Ivanchenko >> 149 } while(i < nbin && x < emax); 145 150 146 if(pVerbose > 1) { << 151 nfilled = i; 147 G4cout << " Probability= " << pProbability << 152 if(fVerbose > 1) { G4cout << " Probability= " << totProbability << G4endl; } 148 << probmax << " emin=" << emin << " << 153 return totProbability; 149 << " E1=" << fE1 << " E2=" << fE2 << G4en << 150 } << 151 return pProbability; << 152 } 154 } 153 155 154 G4double G4VEmissionProbability::SampleEnergy( 156 G4double G4VEmissionProbability::SampleEnergy() 155 { 157 { 156 static const G4double fact = 1.05; << 158 static const G4double fact = 1.1; 157 static const G4double alim = 0.05; << 158 static const G4double blim = 20.; << 159 probmax *= fact; 159 probmax *= fact; >> 160 if(0.0 == fProb[nfilled] && nfilled > 2) { --nfilled; } >> 161 for(size_t i=0; i<=nfilled; ++i) { >> 162 fProb[i] *= fact; >> 163 } >> 164 G4double ekin(0.0), s1(1.0), s2(0.0), ksi(0.0), psi(1.0), z(0.0); 160 165 161 // two regions with flat and exponential maj << 166 if(fVerbose > 1) { 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*probm << 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 << 175 } << 176 p2 = a1/a0; << 177 p1 /= (p1 + p2); << 178 p2 = 1.0 - p1; << 179 } << 180 << 181 if(pVerbose > 1) { << 182 G4cout << "### G4VEmissionProbability::Sam 167 G4cout << "### G4VEmissionProbability::SampleEnergy: " 183 << " Emin= " << emin << " Emax= " << emax 168 << " Emin= " << emin << " Emax= " << emax 184 << "/n E1=" << fE1 << " p1=" << << 169 << " Nf= " << nfilled << G4endl; 185 << " probmax=" << probmax << " P2=" << fP << 170 } >> 171 G4double x0 = emax; >> 172 G4double x1 = fEner[nfilled-1]; >> 173 G4double x2 = fEner[nfilled]; >> 174 G4double y0 = probmax; >> 175 G4double y1 = fProb[nfilled-1]; >> 176 G4double y2 = fProb[nfilled]; >> 177 G4bool islog(false); >> 178 // a condition if a special treatment of falling down >> 179 // distribution is needed >> 180 // >> 181 G4double ee = 0.5*(emax - emin); >> 182 if(nfilled > 5 && y2 > 0.0 && y2 < 0.1*y0 && y1 > 0.0 && x1 - emin < ee) { >> 183 ksi = G4Log(y1/y2)/G4Log(x2/x1); >> 184 x0 = x2*G4Exp(-G4Log(y0/y2)/ksi); >> 185 // general condition to have two sampling area >> 186 if(x0 < x1 && x0 > eprobmax && x0 - emin < ee) { >> 187 // condition when the first area does not exist >> 188 if(x0 <= emin) { >> 189 s1 = 0.0; >> 190 s2 = 1.0; >> 191 y0 *= G4Exp(G4Log(x0/emin)*ksi); >> 192 x0 = emin; >> 193 } else { >> 194 s1 = (x0 - emin)*y0; >> 195 } >> 196 // parameters of the second area >> 197 if(std::abs(1.0 - ksi) < 0.1) { >> 198 z = G4Log(emax/x0); >> 199 if(s1 > 0.0) { s2 = x0*y0*z; } >> 200 islog = true; >> 201 } else { >> 202 psi = 1.0/(1.0 - ksi); >> 203 z = G4Exp(G4Log(emax/x0)*(1. - ksi)) - 1.; >> 204 if(s1 > 0.0) { s2 = std::max(y0*x0*z*psi, 0.0); } >> 205 } >> 206 } >> 207 } >> 208 G4double sum = s1 + s2; >> 209 if(fVerbose > 1) { >> 210 G4cout << " Epeak= " << eprobmax << " e0= " << x0 << " e1= " << x1 << " e2= " << x2 >> 211 << " psi= " << psi << G4endl; >> 212 G4cout << " s1= " << s1 << " s2= " << s2 >> 213 << " y0= " << y0 << " y1= " << y1 << " y2= " << y2 << " z= " << z >> 214 << G4endl; >> 215 for(size_t i=0; i<=nfilled; ++i) { >> 216 G4cout << i << ". E= " << fEner[i] << " P= " << fProb[i] << G4endl; >> 217 } 186 } 218 } 187 219 188 CLHEP::HepRandomEngine* rndm = G4Random::get 220 CLHEP::HepRandomEngine* rndm = G4Random::getTheEngine(); 189 const G4int nmax = 1000; << 221 static const G4int nmax = 100; 190 G4double ekin, gg, gmax; << 222 G4double gr, g; 191 G4int n = 0; << 223 for(size_t i=0; i<nmax; ++i) { 192 do { << 224 G4double q = sum*rndm->flat(); 193 ++n; << 225 if(s2 <= 0.0 || q <= s1) { 194 G4double q = rndm->flat(); << 226 gr = y0; 195 if (p2 == 0.0) { << 227 ekin = emin + (x0 - emin)*q/s1; 196 gmax = probmax; << 228 } else if(islog) { 197 ekin = del*q + emin; << 229 ekin = x0*G4Exp((q - s1)*z/s2); 198 } else if (q <= p1) { << 230 gr = y0*x0/ekin; 199 gmax = probmax; << 200 ekin = del*q/p1 + emin; << 201 } else { 231 } else { 202 ekin = fE1 - G4Log(1.0 - (q - p1)*a1/p2) << 232 ekin = x0*G4Exp(G4Log(1.0 + (q - s1)*z/s2)*psi); 203 x = a0*(ekin - fE1); << 233 gr = y0*G4Exp(G4Log(x0/ekin)*ksi); 204 gmax = fP2; << 205 if(x < blim) { << 206 gmax = probmax*((x > alim) ? G4Exp(-x) : 1.0 << 207 } << 208 } 234 } 209 gg = ComputeProbability(ekin, eCoulomb); << 235 g = ComputeProbability(ekin, eCoulomb); 210 if(pVerbose > 2) { << 236 if(fVerbose > 1) { 211 G4cout << " " << n << 237 G4cout << " " << i 212 << ". prob= " << gg << " probmax= " << << 238 << ". prob= " << g << " probmax= " << gr 213 << " Ekin= " << ekin << G4endl; 239 << " Ekin= " << ekin << G4endl; 214 } 240 } 215 if((gg > gmax || n > nmax) && pVerbose > 1 << 241 if(g > gr && fVerbose > 0) { 216 G4cout << "### G4VEmissionProbability::S << 242 G4cout << "### G4VEmissionProbability::SampleEnergy Warning i= " << i 217 << " A= " << theA << " Eex(MeV)=" << 243 << " prob/probmax= " << g/gr << " rndm= " << q 218 << "\n Warning n= " << n << 244 << "\n prob= " << g << " probmax= " << gr 219 << " prob/gmax=" << gg/gmax << 245 << " z= " << z << " ksi= " << ksi 220 << " prob=" << gg << " gmax=" << gmax < << 221 << "\n Ekin= " << ekin << " Emin= " 246 << "\n Ekin= " << ekin << " Emin= " << emin 222 << " Emax= " << emax << G4endl; << 247 << " Emax= " << emax << " E0= " << x0 << G4endl; >> 248 G4cout << " s1= " << s1 << " s2= " << s2 >> 249 << " y0= " << y0 << " y1= " << y1 << " y2= " << y2 << G4endl; >> 250 for(size_t j=0; j<=nfilled; ++j) { >> 251 G4cout << j << ". E= " << fEner[j] << " P= " << fProb[j] << G4endl; >> 252 } 223 } 253 } 224 } while(gmax*rndm->flat() > gg && n < nmax); << 254 if(gr*rndm->flat() <= g) { break; } 225 G4double enew = FindRecoilExcitation(ekin); << 226 if(pVerbose > 1) { << 227 G4cout << "### SampleEnergy: Efinal= " << 228 << enew << " E=" << ekin << " Eexc=" << << 229 } 255 } 230 return enew; << 256 return ekin; 231 } 257 } 232 258 233 G4double G4VEmissionProbability::FindRecoilExc << 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.*pMa << 241 << 242 fExcRes = mres - pResMass; << 243 << 244 if(pVerbose > 1) { << 245 G4cout << "### FindRecoilExcitation for re << 246 << resZ << " resA= " << resA << 247 << " evaporated Z= " << theZ << " A << 248 << " Ekin= " << e << " Eexc= " << fExcRes << 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)/pMas << 255 } << 256 if(!fFD) { return e; } << 257 << 258 // select final state excitation << 259 auto lManager = pNuclearLevelData->GetLevelM << 260 if(nullptr == lManager) { return e; } << 261 << 262 // levels are not known << 263 if(fExcRes > lManager->MaxLevelEnergy() + pT << 264 << 265 // find level << 266 std::size_t idx = lManager->NearestLevelInde << 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 || p << 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 << 283 } << 284 259