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