[819] | 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 | // |
---|
[1315] | 27 | // $Id: G4PreCompoundEmission.cc,v 1.28 2010/02/25 10:27:36 vnivanch Exp $ |
---|
| 28 | // GEANT4 tag $Name: geant4-09-04-beta-cand-01 $ |
---|
[1055] | 29 | // |
---|
| 30 | // ------------------------------------------------------------------- |
---|
| 31 | // |
---|
| 32 | // GEANT4 Class file |
---|
| 33 | // |
---|
| 34 | // |
---|
| 35 | // File name: G4PreCompoundEmission |
---|
| 36 | // |
---|
| 37 | // Author: V.Lara |
---|
| 38 | // |
---|
| 39 | // Modified: |
---|
[1315] | 40 | // 15.01.2010 J.M.Quesada added protection against unphysical values of parameter an |
---|
| 41 | // 19.01.2010 V.Ivanchenko simplified computation of parameter an, sample cosTheta |
---|
| 42 | // instead of theta; protect all calls to sqrt |
---|
[1055] | 43 | // |
---|
[819] | 44 | |
---|
| 45 | #include "G4PreCompoundEmission.hh" |
---|
| 46 | #include "G4PreCompoundParameters.hh" |
---|
| 47 | |
---|
| 48 | #include "G4PreCompoundEmissionFactory.hh" |
---|
| 49 | #include "G4HETCEmissionFactory.hh" |
---|
| 50 | |
---|
[1315] | 51 | const G4PreCompoundEmission & |
---|
| 52 | G4PreCompoundEmission::operator=(const G4PreCompoundEmission &) |
---|
[819] | 53 | { |
---|
[1315] | 54 | throw G4HadronicException(__FILE__, __LINE__, |
---|
| 55 | "G4PreCompoundEmission::operator= meant to not be accessable"); |
---|
[819] | 56 | return *this; |
---|
| 57 | } |
---|
| 58 | |
---|
| 59 | |
---|
[1055] | 60 | G4bool G4PreCompoundEmission::operator==(const G4PreCompoundEmission &) const |
---|
[819] | 61 | { |
---|
| 62 | return false; |
---|
| 63 | } |
---|
| 64 | |
---|
[1055] | 65 | G4bool G4PreCompoundEmission::operator!=(const G4PreCompoundEmission &) const |
---|
[819] | 66 | { |
---|
| 67 | return true; |
---|
| 68 | } |
---|
| 69 | |
---|
| 70 | G4PreCompoundEmission::G4PreCompoundEmission() |
---|
| 71 | { |
---|
| 72 | theFragmentsFactory = new G4PreCompoundEmissionFactory(); |
---|
| 73 | theFragmentsVector = new G4PreCompoundFragmentVector(theFragmentsFactory->GetFragmentVector()); |
---|
| 74 | } |
---|
| 75 | |
---|
| 76 | G4PreCompoundEmission::~G4PreCompoundEmission() |
---|
| 77 | { |
---|
| 78 | if (theFragmentsFactory) delete theFragmentsFactory; |
---|
| 79 | if (theFragmentsVector) delete theFragmentsVector; |
---|
| 80 | } |
---|
| 81 | |
---|
| 82 | void G4PreCompoundEmission::SetDefaultModel() |
---|
| 83 | { |
---|
| 84 | if (theFragmentsFactory) delete theFragmentsFactory; |
---|
| 85 | theFragmentsFactory = new G4PreCompoundEmissionFactory(); |
---|
| 86 | if (theFragmentsVector) |
---|
| 87 | { |
---|
| 88 | theFragmentsVector->SetVector(theFragmentsFactory->GetFragmentVector()); |
---|
| 89 | } |
---|
| 90 | else |
---|
| 91 | { |
---|
| 92 | theFragmentsVector = new G4PreCompoundFragmentVector(theFragmentsFactory->GetFragmentVector()); |
---|
| 93 | } |
---|
| 94 | theFragmentsVector->ResetStage(); |
---|
| 95 | return; |
---|
| 96 | } |
---|
| 97 | |
---|
| 98 | void G4PreCompoundEmission::SetHETCModel() |
---|
| 99 | { |
---|
| 100 | if (theFragmentsFactory) delete theFragmentsFactory; |
---|
| 101 | theFragmentsFactory = new G4HETCEmissionFactory(); |
---|
| 102 | if (theFragmentsVector) |
---|
| 103 | { |
---|
| 104 | theFragmentsVector->SetVector(theFragmentsFactory->GetFragmentVector()); |
---|
| 105 | } |
---|
| 106 | else |
---|
| 107 | { |
---|
| 108 | theFragmentsVector = new G4PreCompoundFragmentVector(theFragmentsFactory->GetFragmentVector()); |
---|
| 109 | } |
---|
| 110 | theFragmentsVector->ResetStage(); |
---|
| 111 | return; |
---|
| 112 | } |
---|
| 113 | |
---|
| 114 | G4ReactionProduct * G4PreCompoundEmission::PerformEmission(G4Fragment & aFragment) |
---|
| 115 | { |
---|
| 116 | #ifdef debug |
---|
| 117 | G4Fragment InitialState(aFragment); |
---|
| 118 | #endif |
---|
| 119 | // Choose a Fragment for emission |
---|
| 120 | G4VPreCompoundFragment * theFragment = theFragmentsVector->ChooseFragment(); |
---|
| 121 | if (theFragment == 0) |
---|
| 122 | { |
---|
| 123 | G4cerr << "G4PreCompoundEmission::PerformEmission : I couldn't choose a fragment\n" |
---|
| 124 | << "while trying to de-excite\n" |
---|
| 125 | << aFragment << '\n'; |
---|
| 126 | throw G4HadronicException(__FILE__, __LINE__, ""); |
---|
| 127 | } |
---|
| 128 | // Kinetic Energy of emitted fragment |
---|
| 129 | G4double KineticEnergyOfEmittedFragment = theFragment->GetKineticEnergy(aFragment); |
---|
| 130 | |
---|
| 131 | // Calculate the fragment momentum (three vector) |
---|
| 132 | G4ThreeVector momentum = AngularDistribution(theFragment,aFragment,KineticEnergyOfEmittedFragment); |
---|
| 133 | |
---|
| 134 | // Mass of emittef fragment |
---|
| 135 | G4double EmittedMass = theFragment->GetNuclearMass(); |
---|
| 136 | |
---|
| 137 | // Now we can calculate the four momentum |
---|
| 138 | // both options are valid and give the same result but 2nd one is faster |
---|
| 139 | // G4LorentzVector EmittedMomentum(momentum,std::sqrt(momentum.mag2()+EmittedMass*EmittedMass)); |
---|
| 140 | G4LorentzVector EmittedMomentum(momentum,EmittedMass+KineticEnergyOfEmittedFragment); |
---|
| 141 | |
---|
| 142 | // Perform Lorentz boost |
---|
| 143 | EmittedMomentum.boost(aFragment.GetMomentum().boostVector()); |
---|
| 144 | |
---|
| 145 | // Set emitted fragment momentum |
---|
| 146 | theFragment->SetMomentum(EmittedMomentum); |
---|
| 147 | |
---|
| 148 | |
---|
| 149 | // NOW THE RESIDUAL NUCLEUS |
---|
| 150 | // ------------------------ |
---|
| 151 | |
---|
| 152 | // Now the residual nucleus. |
---|
| 153 | // The energy conservation says that |
---|
| 154 | G4double ResidualEcm = |
---|
| 155 | // aFragment.GetGroundStateMass() + aFragment.GetExcitationEnergy() // initial energy in cm |
---|
| 156 | aFragment.GetMomentum().m() |
---|
| 157 | - (EmittedMass+KineticEnergyOfEmittedFragment); |
---|
| 158 | |
---|
| 159 | // Then the four momentum for residual is |
---|
| 160 | G4LorentzVector RestMomentum(-momentum,ResidualEcm); |
---|
| 161 | // This could save a Lorentz boost |
---|
| 162 | // G4LorentzVector RestMomentum2(aFragment.GetMomentum()-EmittedMomentum); |
---|
| 163 | |
---|
| 164 | // Just for test |
---|
| 165 | // Excitation energy |
---|
| 166 | // G4double anU = ResidualEcm - theFragment->GetRestNuclearMass(); |
---|
| 167 | // This is equivalent |
---|
| 168 | // G4double anU = theFragment->GetMaximalKineticEnergy() - KineticEnergyOfEmittedFragment + |
---|
| 169 | // theFragment->GetCoulombBarrier(); |
---|
| 170 | |
---|
| 171 | // check that Excitation energy is >= 0 |
---|
| 172 | G4double anU = RestMomentum.m()-theFragment->GetRestNuclearMass(); |
---|
[1315] | 173 | if (anU < 0.0) { |
---|
| 174 | throw G4HadronicException(__FILE__, __LINE__, |
---|
| 175 | "G4PreCompoundModel::DeExcite: Excitation energy less than 0!"); |
---|
| 176 | } |
---|
| 177 | |
---|
[819] | 178 | // Update nucleus parameters: |
---|
| 179 | // -------------------------- |
---|
[962] | 180 | |
---|
[819] | 181 | // Number of excitons |
---|
| 182 | aFragment.SetNumberOfParticles(aFragment.GetNumberOfParticles()- |
---|
| 183 | static_cast<G4int>(theFragment->GetA())); |
---|
| 184 | // Number of charges |
---|
| 185 | aFragment.SetNumberOfCharged(aFragment.GetNumberOfCharged()- |
---|
| 186 | static_cast<G4int>(theFragment->GetZ())); |
---|
| 187 | |
---|
| 188 | // Atomic number |
---|
| 189 | aFragment.SetA(theFragment->GetRestA()); |
---|
| 190 | |
---|
| 191 | // Charge |
---|
| 192 | aFragment.SetZ(theFragment->GetRestZ()); |
---|
[962] | 193 | |
---|
[819] | 194 | |
---|
| 195 | // Perform Lorentz boosts |
---|
| 196 | RestMomentum.boost(aFragment.GetMomentum().boostVector()); |
---|
| 197 | |
---|
| 198 | // Update nucleus momentum |
---|
| 199 | aFragment.SetMomentum(RestMomentum); |
---|
| 200 | |
---|
| 201 | // Create a G4ReactionProduct |
---|
| 202 | G4ReactionProduct * MyRP = theFragment->GetReactionProduct(); |
---|
| 203 | #ifdef PRECOMPOUND_TEST |
---|
| 204 | MyRP->SetCreatorModel("G4PreCompoundModel"); |
---|
| 205 | #endif |
---|
| 206 | #ifdef debug |
---|
| 207 | CheckConservation(InitialState,aFragment,MyRP); |
---|
| 208 | #endif |
---|
| 209 | return MyRP; |
---|
| 210 | } |
---|
| 211 | |
---|
[1315] | 212 | G4ThreeVector |
---|
| 213 | G4PreCompoundEmission::AngularDistribution(G4VPreCompoundFragment * theFragment, |
---|
| 214 | const G4Fragment& aFragment, |
---|
| 215 | const G4double kinEnergyOfEmittedFrag) const |
---|
[819] | 216 | { |
---|
| 217 | G4double p = aFragment.GetNumberOfParticles(); |
---|
| 218 | G4double h = aFragment.GetNumberOfHoles(); |
---|
| 219 | G4double U = aFragment.GetExcitationEnergy(); |
---|
[1315] | 220 | |
---|
| 221 | G4double ekin = std::max(0.0, kinEnergyOfEmittedFrag); |
---|
[819] | 222 | |
---|
| 223 | // Emission particle separation energy |
---|
| 224 | G4double Bemission = theFragment->GetBindingEnergy(); |
---|
[1315] | 225 | |
---|
[819] | 226 | // Fermi energy |
---|
| 227 | G4double Ef = G4PreCompoundParameters::GetAddress()->GetFermiEnergy(); |
---|
| 228 | |
---|
| 229 | // |
---|
| 230 | // G4EvaporationLevelDensityParameter theLDP; |
---|
| 231 | // G4double g = (6.0/pi2)*aFragment.GetA()* |
---|
[1315] | 232 | |
---|
| 233 | G4double g = (6.0/pi2)*aFragment.GetA() |
---|
| 234 | *G4PreCompoundParameters::GetAddress()->GetLevelDensity(); |
---|
[819] | 235 | |
---|
| 236 | // Average exciton energy relative to bottom of nuclear well |
---|
| 237 | G4double Eav = 2.0*p*(p+1.0)/((p+h)*g); |
---|
| 238 | |
---|
| 239 | // Excitation energy relative to the Fermi Level |
---|
| 240 | G4double Uf = std::max(U - (p - h)*Ef , 0.0); |
---|
| 241 | // G4double Uf = U - KineticEnergyOfEmittedFragment - Bemission; |
---|
| 242 | |
---|
| 243 | G4double w_num = rho(p+1,h,g,Uf,Ef); |
---|
| 244 | G4double w_den = rho(p,h,g,Uf,Ef); |
---|
| 245 | if (w_num > 0.0 && w_den > 0.0) |
---|
| 246 | { |
---|
| 247 | Eav *= (w_num/w_den); |
---|
| 248 | Eav += - Uf/(p+h) + Ef; |
---|
| 249 | } |
---|
| 250 | else |
---|
| 251 | { |
---|
| 252 | Eav = Ef; |
---|
| 253 | } |
---|
| 254 | |
---|
[1196] | 255 | |
---|
[1315] | 256 | // VI + JMQ 19/01/2010 update computation of the parameter an |
---|
| 257 | // |
---|
| 258 | G4double an = 0.0; |
---|
| 259 | G4double Eeff = ekin + Bemission + Ef; |
---|
| 260 | if(ekin > DBL_MIN && Eeff > DBL_MIN) { |
---|
| 261 | |
---|
| 262 | G4double zeta = std::max(1.0,9.3/std::sqrt(ekin/MeV)); |
---|
| 263 | |
---|
| 264 | an = 3.0*std::sqrt((ProjEnergy+Ef)*Eeff)/(zeta*Eav); |
---|
| 265 | |
---|
| 266 | G4int ne = aFragment.GetNumberOfExcitons() - 1; |
---|
| 267 | if ( ne > 1 ) { an /= (G4double)ne; } |
---|
[819] | 268 | |
---|
[1315] | 269 | // protection of exponent |
---|
| 270 | if ( an > 10. ) { an = 10.; } |
---|
[1196] | 271 | } |
---|
| 272 | |
---|
[1315] | 273 | // sample cosine of theta and not theta as in old versions |
---|
| 274 | G4double random = G4UniformRand(); |
---|
| 275 | G4double cost; |
---|
| 276 | |
---|
| 277 | if(an < 0.1) { cost = 1. - 2*random; } |
---|
| 278 | else { |
---|
| 279 | G4double exp2an = std::exp(-2*an); |
---|
| 280 | cost = 1. + std::log(1-random*(1-exp2an))/an; |
---|
| 281 | if(cost > 1.) { cost = 1.; } |
---|
| 282 | else if(cost < -1.) {cost = -1.; } |
---|
| 283 | } |
---|
| 284 | |
---|
[819] | 285 | G4double phi = twopi*G4UniformRand(); |
---|
| 286 | |
---|
| 287 | // Calculate the momentum magnitude of emitted fragment |
---|
[1315] | 288 | G4double pmag = std::sqrt(ekin*(ekin + 2.0*theFragment->GetNuclearMass())); |
---|
[819] | 289 | |
---|
[1315] | 290 | G4double sint = std::sqrt((1.0-cost)*(1.0+cost)); |
---|
[819] | 291 | |
---|
[1315] | 292 | G4ThreeVector momentum(pmag*std::cos(phi)*sint,pmag*std::sin(phi)*sint,pmag*cost); |
---|
[819] | 293 | // theta is the angle wrt the incident direction |
---|
| 294 | momentum.rotateUz(theIncidentDirection); |
---|
| 295 | |
---|
| 296 | return momentum; |
---|
| 297 | } |
---|
| 298 | |
---|
| 299 | G4double G4PreCompoundEmission::rho(const G4double p, const G4double h, const G4double g, |
---|
| 300 | const G4double E, const G4double Ef) const |
---|
[1315] | 301 | { |
---|
| 302 | // 25.02.2010 V.Ivanchenko added more protections |
---|
| 303 | G4double Aph = (p*p + h*h + p - 3.0*h)/(4.0*g); |
---|
| 304 | // G4double alpha = (p*p + h*h)/(2.0*g); |
---|
[962] | 305 | |
---|
[1315] | 306 | if ( E - Aph < 0.0) { return 0.0; } |
---|
[962] | 307 | |
---|
[1196] | 308 | G4double logConst = (p+h)*std::log(g) - logfactorial(p+h-1) - logfactorial(p) - logfactorial(h); |
---|
[962] | 309 | |
---|
[1315] | 310 | // initialise values using j=0 |
---|
[962] | 311 | |
---|
| 312 | G4double t1=1; |
---|
| 313 | G4double t2=1; |
---|
[1315] | 314 | G4double logt3=(p+h-1) * std::log(E-Aph) + logConst; |
---|
| 315 | const G4double logmax = 200.; |
---|
| 316 | if(logt3 > logmax) { logt3 = logmax; } |
---|
| 317 | G4double tot = std::exp( logt3 ); |
---|
[962] | 318 | |
---|
[1315] | 319 | // and now sum rest of terms |
---|
| 320 | // 25.02.2010 V.Ivanchenko change while to for loop and cleanup |
---|
| 321 | G4double Eeff = E - Aph; |
---|
| 322 | for(G4int j=1; j<=h; ++j) |
---|
[962] | 323 | { |
---|
[1315] | 324 | Eeff -= Ef; |
---|
| 325 | if(Eeff < 0.0) { break; } |
---|
| 326 | t1 *= -1.; |
---|
| 327 | t2 *= (G4double)(h+1-j)/(G4double)j; |
---|
| 328 | logt3 = (p+h-1) * std::log( Eeff) + logConst; |
---|
| 329 | if(logt3 > logmax) { logt3 = logmax; } |
---|
| 330 | tot += t1*t2*std::exp(logt3); |
---|
[962] | 331 | } |
---|
| 332 | |
---|
| 333 | return tot; |
---|
| 334 | } |
---|
| 335 | |
---|
| 336 | G4double G4PreCompoundEmission::factorial(G4double a) const |
---|
| 337 | { |
---|
[819] | 338 | // Values of factorial function from 0 to 60 |
---|
| 339 | const G4int factablesize = 61; |
---|
| 340 | static const G4double fact[factablesize] = |
---|
| 341 | { |
---|
| 342 | 1.0, // 0! |
---|
| 343 | 1.0, // 1! |
---|
| 344 | 2.0, // 2! |
---|
| 345 | 6.0, // 3! |
---|
| 346 | 24.0, // 4! |
---|
| 347 | 120.0, // 5! |
---|
| 348 | 720.0, // 6! |
---|
| 349 | 5040.0, // 7! |
---|
| 350 | 40320.0, // 8! |
---|
| 351 | 362880.0, // 9! |
---|
| 352 | 3628800.0, // 10! |
---|
| 353 | 39916800.0, // 11! |
---|
| 354 | 479001600.0, // 12! |
---|
| 355 | 6227020800.0, // 13! |
---|
| 356 | 87178291200.0, // 14! |
---|
| 357 | 1307674368000.0, // 15! |
---|
| 358 | 20922789888000.0, // 16! |
---|
| 359 | 355687428096000.0, // 17! |
---|
| 360 | 6402373705728000.0, // 18! |
---|
| 361 | 121645100408832000.0, // 19! |
---|
| 362 | 2432902008176640000.0, // 20! |
---|
| 363 | 51090942171709440000.0, // 21! |
---|
| 364 | 1124000727777607680000.0, // 22! |
---|
| 365 | 25852016738884976640000.0, // 23! |
---|
| 366 | 620448401733239439360000.0, // 24! |
---|
| 367 | 15511210043330985984000000.0, // 25! |
---|
| 368 | 403291461126605635584000000.0, // 26! |
---|
| 369 | 10888869450418352160768000000.0, // 27! |
---|
| 370 | 304888344611713860501504000000.0, // 28! |
---|
| 371 | 8841761993739701954543616000000.0, // 29! |
---|
| 372 | 265252859812191058636308480000000.0, // 30! |
---|
| 373 | 8222838654177922817725562880000000.0, // 31! |
---|
| 374 | 263130836933693530167218012160000000.0, // 32! |
---|
| 375 | 8683317618811886495518194401280000000.0, // 33! |
---|
| 376 | 295232799039604140847618609643520000000.0, // 34! |
---|
| 377 | 10333147966386144929666651337523200000000.0, // 35! |
---|
| 378 | 371993326789901217467999448150835200000000.0, // 36! |
---|
| 379 | 13763753091226345046315979581580902400000000.0, // 37! |
---|
| 380 | 523022617466601111760007224100074291200000000.0, // 38! |
---|
| 381 | 20397882081197443358640281739902897356800000000.0, // 39! |
---|
| 382 | 815915283247897734345611269596115894272000000000.0, // 40! |
---|
| 383 | 33452526613163807108170062053440751665152000000000.0, // 41! |
---|
| 384 | 1405006117752879898543142606244511569936384000000000.0, // 42! |
---|
| 385 | 60415263063373835637355132068513997507264512000000000.0, // 43! |
---|
| 386 | 2658271574788448768043625811014615890319638528000000000.0, // 44! |
---|
| 387 | 119622220865480194561963161495657715064383733760000000000.0, // 45! |
---|
| 388 | 5502622159812088949850305428800254892961651752960000000000.0, // 46! |
---|
| 389 | 258623241511168180642964355153611979969197632389120000000000.0, // 47! |
---|
| 390 | 12413915592536072670862289047373375038521486354677760000000000.0, // 48! |
---|
| 391 | 608281864034267560872252163321295376887552831379210240000000000.0, // 49! |
---|
| 392 | 30414093201713378043612608166064768844377641568960512000000000000.0, // 50! |
---|
| 393 | 1551118753287382280224243016469303211063259720016986112000000000000.0, // 51! |
---|
| 394 | 80658175170943878571660636856403766975289505440883277824000000000000.0, // 52! |
---|
| 395 | 4274883284060025564298013753389399649690343788366813724672000000000000.0, // 53! |
---|
| 396 | 230843697339241380472092742683027581083278564571807941132288000000000000.0, // 54! |
---|
| 397 | 12696403353658275925965100847566516959580321051449436762275840000000000000.0, // 55! |
---|
| 398 | 710998587804863451854045647463724949736497978881168458687447040000000000000.0, // 56! |
---|
| 399 | 40526919504877216755680601905432322134980384796226602145184481280000000000000.0, // 57! |
---|
| 400 | 2350561331282878571829474910515074683828862318181142924420699914240000000000000.0, // 58! |
---|
| 401 | 138683118545689835737939019720389406345902876772687432540821294940160000000000000.0, // 59! |
---|
| 402 | 8320987112741390144276341183223364380754172606361245952449277696409600000000000000.0 // 60! |
---|
| 403 | }; |
---|
| 404 | // fact[0] = 1; |
---|
| 405 | // for (G4int n = 1; n < 21; n++) { |
---|
| 406 | // fact[n] = fact[n-1]*static_cast<G4double>(n); |
---|
| 407 | // } |
---|
[962] | 408 | G4double result(0.0); |
---|
| 409 | G4int ia = static_cast<G4int>(a); |
---|
| 410 | if (ia < factablesize) |
---|
[819] | 411 | { |
---|
[962] | 412 | result = fact[ia]; |
---|
[819] | 413 | } |
---|
| 414 | else |
---|
| 415 | { |
---|
[962] | 416 | result = fact[factablesize-1]; |
---|
| 417 | for (G4int n = factablesize; n < ia+1; ++n) |
---|
[819] | 418 | { |
---|
[962] | 419 | result *= static_cast<G4double>(n); |
---|
[819] | 420 | } |
---|
| 421 | } |
---|
| 422 | |
---|
[962] | 423 | return result; |
---|
[819] | 424 | } |
---|
[1196] | 425 | G4double G4PreCompoundEmission::logfactorial(G4double a) const |
---|
| 426 | { |
---|
| 427 | // Values of logs of factorial function from 0 to 60 |
---|
[819] | 428 | |
---|
[1196] | 429 | G4double result(0.0); |
---|
| 430 | const G4int factablesize = 61; |
---|
| 431 | const G4double halfLn2pi = 0.918938533; // 0.5 log(2* pi) |
---|
| 432 | static G4double logfact[factablesize]; |
---|
| 433 | static bool needinit=true; |
---|
| 434 | |
---|
| 435 | if (needinit) |
---|
| 436 | { |
---|
| 437 | needinit=false; |
---|
| 438 | for ( G4int n=0; n < factablesize; ++n) |
---|
| 439 | { |
---|
| 440 | logfact[n]=std::log(factorial(n)); |
---|
| 441 | } |
---|
| 442 | } |
---|
[819] | 443 | |
---|
[1196] | 444 | G4int ia = static_cast<G4int>(a); |
---|
| 445 | if (ia < factablesize) |
---|
| 446 | { |
---|
| 447 | result = logfact[ia]; |
---|
| 448 | } else { |
---|
| 449 | result = (ia+0.5)*std::log(G4double(ia)) - ia + halfLn2pi; |
---|
| 450 | } |
---|
| 451 | |
---|
| 452 | return result; |
---|
| 453 | } |
---|
| 454 | |
---|
[819] | 455 | #ifdef debug |
---|
| 456 | void G4PreCompoundEmission::CheckConservation(const G4Fragment & theInitialState, |
---|
| 457 | const G4Fragment & theResidual, |
---|
| 458 | G4ReactionProduct * theEmitted) const |
---|
| 459 | { |
---|
| 460 | G4double ProductsEnergy = theEmitted->GetTotalEnergy() + theResidual.GetMomentum().e(); |
---|
| 461 | G4ThreeVector ProductsMomentum(theEmitted->GetMomentum()+theResidual.GetMomentum().vect()); |
---|
| 462 | G4int ProductsA = theEmitted->GetDefinition()->GetBaryonNumber() + theResidual.GetA(); |
---|
| 463 | G4int ProductsZ = theEmitted->GetDefinition()->GetPDGCharge() + theResidual.GetZ(); |
---|
| 464 | |
---|
| 465 | if (ProductsA != theInitialState.GetA()) { |
---|
| 466 | G4cout << "!!!!!!!!!! Baryonic Number Conservation Violation !!!!!!!!!!" << G4endl; |
---|
| 467 | G4cout << "G4PreCompoundEmission.cc: Barionic Number Conservation" |
---|
| 468 | << G4endl; |
---|
| 469 | G4cout << "Initial A = " << theInitialState.GetA() |
---|
| 470 | << " Fragments A = " << ProductsA << " Diference --> " |
---|
| 471 | << theInitialState.GetA() - ProductsA << G4endl; |
---|
| 472 | } |
---|
| 473 | if (ProductsZ != theInitialState.GetZ()) { |
---|
| 474 | G4cout << "!!!!!!!!!! Charge Conservation Violation !!!!!!!!!!" << G4endl; |
---|
| 475 | G4cout << "G4PreCompoundEmission.cc: Charge Conservation test" |
---|
| 476 | << G4endl; |
---|
| 477 | G4cout << "Initial Z = " << theInitialState.GetZ() |
---|
| 478 | << " Fragments Z = " << ProductsZ << " Diference --> " |
---|
| 479 | << theInitialState.GetZ() - ProductsZ << G4endl; |
---|
| 480 | } |
---|
| 481 | if (std::abs(ProductsEnergy-theInitialState.GetMomentum().e()) > 10.0*eV) { |
---|
| 482 | G4cout << "!!!!!!!!!! Energy Conservation Violation !!!!!!!!!!" << G4endl; |
---|
| 483 | G4cout << "G4PreCompoundEmission.cc: Energy Conservation test" |
---|
| 484 | << G4endl; |
---|
| 485 | G4cout << "Initial E = " << theInitialState.GetMomentum().e()/MeV << " MeV" |
---|
| 486 | << " Fragments E = " << ProductsEnergy/MeV << " MeV Diference --> " |
---|
| 487 | << (theInitialState.GetMomentum().e() - ProductsEnergy)/MeV << " MeV" << G4endl; |
---|
| 488 | } |
---|
| 489 | if (std::abs(ProductsMomentum.x()-theInitialState.GetMomentum().x()) > 10.0*eV || |
---|
| 490 | std::abs(ProductsMomentum.y()-theInitialState.GetMomentum().y()) > 10.0*eV || |
---|
| 491 | std::abs(ProductsMomentum.z()-theInitialState.GetMomentum().z()) > 10.0*eV) { |
---|
| 492 | G4cout << "!!!!!!!!!! Momentum Conservation Violation !!!!!!!!!!" << G4endl; |
---|
| 493 | G4cout << "G4PreCompoundEmission.cc: Momentum Conservation test" |
---|
| 494 | << G4endl; |
---|
| 495 | G4cout << "Initial P = " << theInitialState.GetMomentum().vect() << " MeV" |
---|
| 496 | << " Fragments P = " << ProductsMomentum << " MeV Diference --> " |
---|
| 497 | << theInitialState.GetMomentum().vect() - ProductsMomentum << " MeV" << G4endl; |
---|
| 498 | } |
---|
| 499 | return; |
---|
| 500 | } |
---|
| 501 | |
---|
| 502 | #endif |
---|