source: trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundEmission.cc @ 847

Last change on this file since 847 was 819, checked in by garnier, 16 years ago

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