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 | // $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 | // |
---|
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 | |
---|
55 | G4bool G4PreCompoundEmission::operator==(const G4PreCompoundEmission &) const |
---|
56 | { |
---|
57 | return false; |
---|
58 | } |
---|
59 | |
---|
60 | G4bool G4PreCompoundEmission::operator!=(const G4PreCompoundEmission &) const |
---|
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 | // -------------------------- |
---|
174 | |
---|
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()); |
---|
187 | |
---|
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 | { |
---|
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 | { |
---|
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 | // } |
---|
400 | G4double result(0.0); |
---|
401 | G4int ia = static_cast<G4int>(a); |
---|
402 | if (ia < factablesize) |
---|
403 | { |
---|
404 | result = fact[ia]; |
---|
405 | } |
---|
406 | else |
---|
407 | { |
---|
408 | result = fact[factablesize-1]; |
---|
409 | for (G4int n = factablesize; n < ia+1; ++n) |
---|
410 | { |
---|
411 | result *= static_cast<G4double>(n); |
---|
412 | } |
---|
413 | } |
---|
414 | |
---|
415 | return result; |
---|
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 |
---|