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

Last change on this file since 1228 was 1228, checked in by garnier, 14 years ago

update geant4.9.3 tag
File size: 18.0 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// $Id: G4PreCompoundEmission.cc,v 1.22 2009/11/13 17:40:14 gcosmo Exp $
28// GEANT4 tag $Name: geant4-09-03 $
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
48const G4PreCompoundEmission & G4PreCompoundEmission::operator=(const G4PreCompoundEmission &)
49{
50  throw G4HadronicException(__FILE__, __LINE__, "G4PreCompoundEmission::operator= meant to not be accessable");
51  return *this;
52}
53
54
55G4bool G4PreCompoundEmission::operator==(const G4PreCompoundEmission &) const
56{
57  return false;
58}
59
60G4bool G4PreCompoundEmission::operator!=(const G4PreCompoundEmission &) const
61{
62  return true;
63}
64
65G4PreCompoundEmission::G4PreCompoundEmission()
66{
67  theFragmentsFactory = new G4PreCompoundEmissionFactory();
68  theFragmentsVector = new G4PreCompoundFragmentVector(theFragmentsFactory->GetFragmentVector());
69}
70
71G4PreCompoundEmission::~G4PreCompoundEmission()
72{
73  if (theFragmentsFactory) delete theFragmentsFactory;
74  if (theFragmentsVector) delete theFragmentsVector;
75}
76
77void 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
93void 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
109G4ReactionProduct * 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
207G4ThreeVector 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
253  G4double zeta = std::max(1.0,9.3/std::sqrt(KineticEnergyOfEmittedFragment/MeV));
254       
255  G4double an = 3.0*std::sqrt((ProjEnergy+Ef)*(KineticEnergyOfEmittedFragment+Bemission+Ef));
256  if (aFragment.GetNumberOfExcitons() == 1)
257    {
258      an /= (zeta*2.0*aFragment.GetNumberOfExcitons()*Eav);
259    }
260  else
261    {
262      an /= (zeta*(aFragment.GetNumberOfExcitons()-1.0)*Eav);
263    }
264                       
265                       
266//    G4double expan=std::exp(an);
267//    thetaold = std::acos(std::log(expan-random*(expan-1.0/expan))/an);
268//      exp(an) becomes large for large an
269//      1/exp(an) becomes large for large negative an
270//       take the large value out of the log depending on the sign of an to avoid FP exceptions
271 
272  G4double random=G4UniformRand();
273  G4double exp2an=0;
274  G4double theta;
275  if (an > 0.) {
276    if (an < 25.) { exp2an = std::exp(-2*an); } // we subtract from 1, exp(-50)~1e-21, so will not
277                                           //  change numerical result
278    theta = std::acos(1+ std::log(1-random*(1-exp2an))/an);
279  } else if ( an < 0.) {
280    if ( an > -25.) { exp2an = std::exp(2*an); } // similar to above, except we compare to rndm*1
281    theta = std::acos(std::log(exp2an-random*(exp2an-1))/an - 1.);
282  } else {   // an==0 now.     
283    theta=std::acos(1.-2*random);
284  }
285 
286  if (std::abs(an) < 50 )
287  {
288     
289  }
290
291  G4double phi = twopi*G4UniformRand();
292 
293  // Calculate the momentum magnitude of emitted fragment       
294  G4double EmittedMass = theFragment->GetNuclearMass();
295  G4double pmag = std::sqrt(KineticEnergyOfEmittedFragment*(KineticEnergyOfEmittedFragment+2.0*EmittedMass));
296 
297 
298  G4double sinTheta = std::sin(theta);
299  //  G4double cosTheta = std::sqrt(1.0-sinTheta*sinTheta);
300  G4double cosTheta = std::cos(theta);
301
302  G4ThreeVector momentum(pmag*std::cos(phi)*sinTheta,pmag*std::sin(phi)*sinTheta,pmag*cosTheta);
303  // theta is the angle wrt the incident direction
304  momentum.rotateUz(theIncidentDirection);
305
306  return momentum;
307}
308
309G4double G4PreCompoundEmission::rho(const G4double p, const G4double h, const G4double g, 
310                                    const G4double E, const G4double Ef) const
311{
312       
313  G4double Aph = (p*p + h*h + p - 3.0*h)/(4.0*g);
314  G4double alpha = (p*p+h*h)/(2.0*g);
315 
316  if ( (E-alpha) < 0 ) return 0;
317 
318  G4double logConst =  (p+h)*std::log(g) - logfactorial(p+h-1) - logfactorial(p) - logfactorial(h);
319
320// initialise values using j=0
321
322  G4double t1=1;
323  G4double t2=1;
324  G4double logt3=(p+h-1) * std::log(E-Aph);
325  G4double tot = std::exp( logt3 + logConst );
326
327// and now sum rest of terms
328  G4int j(1); 
329  while ( (j <= h) && ((E - alpha - j*Ef) > 0.0) ) 
330    {
331          t1 *= -1.;
332          t2 *= (h+1-j)/j;
333          logt3 = (p+h-1) * std::log( E - j*Ef - Aph) + logConst;
334          G4double t3 = std::exp(logt3);
335          tot += t1*t2*t3;
336          j++;
337    }
338       
339  return tot;
340}
341
342
343
344G4double G4PreCompoundEmission::factorial(G4double a) const
345{
346  // Values of factorial function from 0 to 60
347  const G4int factablesize = 61;
348  static const G4double fact[factablesize] = 
349    {
350      1.0, // 0!
351      1.0, // 1!
352      2.0, // 2!
353      6.0, // 3!
354      24.0, // 4!
355      120.0, // 5!
356      720.0, // 6!
357      5040.0, // 7!
358      40320.0, // 8!
359      362880.0, // 9!
360      3628800.0, // 10!
361      39916800.0, // 11!
362      479001600.0, // 12!
363      6227020800.0, // 13!
364      87178291200.0, // 14!
365      1307674368000.0, // 15!
366      20922789888000.0, // 16!
367      355687428096000.0, // 17!
368      6402373705728000.0, // 18!
369      121645100408832000.0, // 19!
370      2432902008176640000.0, // 20!
371      51090942171709440000.0, // 21!
372      1124000727777607680000.0, // 22!
373      25852016738884976640000.0, // 23!
374      620448401733239439360000.0, // 24!
375      15511210043330985984000000.0, // 25!
376      403291461126605635584000000.0, // 26!
377      10888869450418352160768000000.0, // 27!
378      304888344611713860501504000000.0, // 28!
379      8841761993739701954543616000000.0, // 29!
380      265252859812191058636308480000000.0, // 30!
381      8222838654177922817725562880000000.0, // 31!
382      263130836933693530167218012160000000.0, // 32!
383      8683317618811886495518194401280000000.0, // 33!
384      295232799039604140847618609643520000000.0, // 34!
385      10333147966386144929666651337523200000000.0, // 35!
386      371993326789901217467999448150835200000000.0, // 36!
387      13763753091226345046315979581580902400000000.0, // 37!
388      523022617466601111760007224100074291200000000.0, // 38!
389      20397882081197443358640281739902897356800000000.0, // 39!
390      815915283247897734345611269596115894272000000000.0, // 40!
391      33452526613163807108170062053440751665152000000000.0, // 41!
392      1405006117752879898543142606244511569936384000000000.0, // 42!
393      60415263063373835637355132068513997507264512000000000.0, // 43!
394      2658271574788448768043625811014615890319638528000000000.0, // 44!
395      119622220865480194561963161495657715064383733760000000000.0, // 45!
396      5502622159812088949850305428800254892961651752960000000000.0, // 46!
397      258623241511168180642964355153611979969197632389120000000000.0, // 47!
398      12413915592536072670862289047373375038521486354677760000000000.0, // 48!
399      608281864034267560872252163321295376887552831379210240000000000.0, // 49!
400      30414093201713378043612608166064768844377641568960512000000000000.0, // 50!
401      1551118753287382280224243016469303211063259720016986112000000000000.0, // 51!
402      80658175170943878571660636856403766975289505440883277824000000000000.0, // 52!
403      4274883284060025564298013753389399649690343788366813724672000000000000.0, // 53!
404      230843697339241380472092742683027581083278564571807941132288000000000000.0, // 54!
405      12696403353658275925965100847566516959580321051449436762275840000000000000.0, // 55!
406      710998587804863451854045647463724949736497978881168458687447040000000000000.0, // 56!
407      40526919504877216755680601905432322134980384796226602145184481280000000000000.0, // 57!
408      2350561331282878571829474910515074683828862318181142924420699914240000000000000.0, // 58!
409      138683118545689835737939019720389406345902876772687432540821294940160000000000000.0, // 59!
410      8320987112741390144276341183223364380754172606361245952449277696409600000000000000.0  // 60!
411    };
412  //    fact[0] = 1;
413  //    for (G4int n = 1; n < 21; n++) {
414  //      fact[n] = fact[n-1]*static_cast<G4double>(n);
415  //    }
416  G4double result(0.0);
417  G4int ia = static_cast<G4int>(a);
418  if (ia < factablesize) 
419    {
420      result = fact[ia];
421    }
422  else
423    {
424      result = fact[factablesize-1];
425      for (G4int n = factablesize; n < ia+1; ++n)
426        {
427          result *= static_cast<G4double>(n);
428        }
429    }
430   
431    return result;
432}
433G4double G4PreCompoundEmission::logfactorial(G4double a) const
434{
435  // Values of logs of factorial function from 0 to 60
436
437  G4double result(0.0);
438  const G4int factablesize = 61;
439  const G4double halfLn2pi = 0.918938533;      // 0.5 log(2* pi)
440  static G4double logfact[factablesize];
441  static bool needinit=true;
442 
443  if (needinit)
444  {
445      needinit=false;
446      for ( G4int n=0; n < factablesize; ++n)
447      {
448         logfact[n]=std::log(factorial(n));
449      }
450  }
451
452  G4int ia = static_cast<G4int>(a);
453  if (ia < factablesize)
454  {
455      result = logfact[ia];
456  } else {
457      result = (ia+0.5)*std::log(G4double(ia)) - ia + halfLn2pi;
458  }
459   
460  return result;
461}
462
463#ifdef debug
464void G4PreCompoundEmission::CheckConservation(const G4Fragment & theInitialState,
465                                              const G4Fragment & theResidual,
466                                              G4ReactionProduct * theEmitted) const
467{
468  G4double ProductsEnergy = theEmitted->GetTotalEnergy() + theResidual.GetMomentum().e();
469  G4ThreeVector ProductsMomentum(theEmitted->GetMomentum()+theResidual.GetMomentum().vect());
470  G4int ProductsA = theEmitted->GetDefinition()->GetBaryonNumber() + theResidual.GetA();
471  G4int ProductsZ = theEmitted->GetDefinition()->GetPDGCharge() + theResidual.GetZ();
472
473  if (ProductsA != theInitialState.GetA()) {
474    G4cout << "!!!!!!!!!! Baryonic Number Conservation Violation !!!!!!!!!!" << G4endl;
475    G4cout << "G4PreCompoundEmission.cc: Barionic Number Conservation"
476           << G4endl; 
477    G4cout << "Initial A = " << theInitialState.GetA() 
478           << "   Fragments A = " << ProductsA << "   Diference --> " 
479           << theInitialState.GetA() - ProductsA << G4endl;
480  }
481  if (ProductsZ != theInitialState.GetZ()) {
482    G4cout << "!!!!!!!!!! Charge Conservation Violation !!!!!!!!!!" << G4endl;
483    G4cout << "G4PreCompoundEmission.cc: Charge Conservation test"
484           << G4endl; 
485    G4cout << "Initial Z = " << theInitialState.GetZ() 
486           << "   Fragments Z = " << ProductsZ << "   Diference --> " 
487           << theInitialState.GetZ() - ProductsZ << G4endl;
488  }
489  if (std::abs(ProductsEnergy-theInitialState.GetMomentum().e()) > 10.0*eV) {
490    G4cout << "!!!!!!!!!! Energy Conservation Violation !!!!!!!!!!" << G4endl;
491    G4cout << "G4PreCompoundEmission.cc: Energy Conservation test" 
492           << G4endl; 
493    G4cout << "Initial E = " << theInitialState.GetMomentum().e()/MeV << " MeV"
494           << "   Fragments E = " << ProductsEnergy/MeV  << " MeV   Diference --> " 
495           << (theInitialState.GetMomentum().e() - ProductsEnergy)/MeV << " MeV" << G4endl;
496  } 
497  if (std::abs(ProductsMomentum.x()-theInitialState.GetMomentum().x()) > 10.0*eV || 
498      std::abs(ProductsMomentum.y()-theInitialState.GetMomentum().y()) > 10.0*eV ||
499      std::abs(ProductsMomentum.z()-theInitialState.GetMomentum().z()) > 10.0*eV) {
500    G4cout << "!!!!!!!!!! Momentum Conservation Violation !!!!!!!!!!" << G4endl;
501    G4cout << "G4PreCompoundEmission.cc: Momentum Conservation test" 
502           << G4endl; 
503    G4cout << "Initial P = " << theInitialState.GetMomentum().vect() << " MeV"
504           << "   Fragments P = " << ProductsMomentum  << " MeV   Diference --> " 
505           << theInitialState.GetMomentum().vect() - ProductsMomentum << " MeV" << G4endl;
506  }
507  return;
508}
509
510#endif
Note: See TracBrowser for help on using the repository browser.